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Executive  Summary 

This  report  was  prepared  by  MDH  Engineered  Solutions  Corp.  (MDH)  to  document 
Phase  II  of  a  project  for  Alberta  Environment  (AE)  to  evaluate  computer  codes  for 
predicting  the  fate  and  transport  of  salt  (sodium  chloride)  to  facilitate  risk-based 
corrective  action. 

Phase  I  of  the  project  recommended  five  modelling  codes  for  detailed  analysis  in 
Phase  II.  The  five  codes  were  chosen  from  a  list  of  almost  250  codes  that  were  ranked 
on  their  strengths,  weaknesses,  and  applicability  to  the  AE  salt  release  scenarios.  The 
five  codes  recommended  for  detailed  analysis  in  Phase  II  were: 

1)  VS2DTI  (USGS); 

2)  UNSATCHEM  (USSL); 

3)  HYDRUS-2D  (USSL-IGWMC); 

4)  SEVIEW  (ESCI);  and, 

5)  CHEMFLOW2000  (OSU). 

The  principal  objective  of  Phase  II  was  to  document  the  strengths  and  weaknesses  of 
each  code  for  the  specific  task  of  simulating  the  transport  of  salt  (sodium  chloride)  at 
relatively  low  concentrations,  in  soil  and  groundwater  in  climatic  regimes  characteristic  of 
the  Alberta  foothills  and  prairies.  Relatively  low  concentrations,  in  this  context,  is 
interpreted  as  concentrations  below  which  density-dependent  behaviour  can  be  ignored. 
This  report  describes  the  performance  of  the  five  codes  for  a  variety  of  generic  scenarios 
with  varying  degrees  of  complexity.  The  generic  scenarios  were  derived  from  a  review 
of  typical  case  histories  and  are  used  to  illustrate  and  compare  the  performance  of  each 
of  the  five  selected  computer  codes.  The  effects  of  density-dependent  flow  and  flow 
through  discretely  fractured  media  were  not  evaluated. 

The  results  of  the  detailed  analysis  indicated  four  of  the  codes  (CHEMFLO,  HYDRUS, 
UNSATCHEM,  and  VS2DTI)  can  be  readily  applied  to  1D  Tier  2  analysis  for  the 
fate/transport  of  salt.  All  four  codes  produced  comparable  consistent  results  for  the 
Tier  2  generic  scenarios. 
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SEVIEW  (AT  123D,  BIOSCREEN  and  SESOIL)  was  the  only  code  tested  that  was 
designed  to  predict  groundwater  recharge  from  seasonal  climatic  data.  SESOIL  is  not 
recommended  because  of  the  difficulties  encountered  in  obtaining  reasonable  estimates 
of  groundwater  recharge  through  tills  for  the  arid  western  Prairie  environment. 
Nevertheless,  AT123D  can  be  used  to  carry  out  1D  Tier  2  analysis  but  without  SESOIL 
has  no  advantages  over  the  other  codes  and  is  limited  to  homogeneous  (or  equivalent 
homogeneous)  materials. 

For  the  more  complex  Tier  3  analysis  only  VS2DTI  and  HYDRUS-2D  were  capable  of 
handling  2D  planar  and  axisymmetric  problems.  Both  codes  performed  consistently  on 
the  2D  generic  scenarios  but  VS2DTI  proved  more  flexible  in  use. 

For  cation  exchange  problems  only  UNSATCHEM  and  VS2DTI  provided  a  means  of 
calculation.  The  two  codes  use  different  definitions  for  selectivity  coefficients  and  are 
difficult  to  compare  head-to-head.  UNSATCHEM  is  the  most  sophisticated  reactive 
transport  model  and  readily  predicts  soil  quality  parameters  such  as  sodium  absorption 
ratio  (SAR)  and  electrical  conductivity  (EC). 

The  detailed  comparison  resulted  in  the  following  recommendations: 

1)  The  revised  version  of  CHEMFLO  to  be  released  after  March  2003  is 
recommended  as  first  choice  for  Tier  2  analysis. 

2)  UNSATCHEM  is  recommended  for  use  in  Tier  3  assessment  where 
simulation  of  cation  exchange  processes  and  estimates  of  SAR  are 
requirements. 

3)  VS2DTI  is  recommended  for  application  to  both  1D  column  (Tier  2)  and  2D 
planar  and  axisymmetric  (Tier  3)  scenarios. 

4)  SEVIEW  is  not  recommended  as  a  suitable  code  for  screening  of  salt 
contaminated  sites  in  Alberta. 

5)  HYDRUS-2D  is  judged  to  be  less  effective  than  VS2DI  as  a  tool  for  Tier  3 


analysis. 
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1.0  INTRODUCTION 

This  report  was  prepared  by  MDH  Engineered  Solutions  Corp.  (MDH)  to  document 
Phase  II  of  a  project  for  Alberta  Environment  (AE)  to  evaluate  computer  codes  for 
predicting  the  fate  and  transport  of  salt  (sodium  chloride)  to  facilitate  risk-based 
corrective  action. 

During  Phase  I  of  the  project  (MDH,  2002),  a  two-level  ranking  procedure  was  used  to 
compare  almost  250  (247)  codes.  The  first  stage  of  the  process  was  based  on  review  of 
readily  available  documentation,  some  discussion  with  model  developers,  and  in-house 
expertise  and  experience  of  MDH  personnel.  This  pre-screening  stage  eliminated  over 
200  (204)  codes,  leaving  43  for  detailed  analysis.  The  remaining  codes  were  ranked 
based  on  17  objective  criteria.  The  top  five  codes  from  the  ranking  matrix  were  reviewed 
in  detail  for  their  strengths,  weaknesses,  and  applicability  to  the  AE  salt  release 
scenarios.  The  five  codes  recommended  for  detailed  analysis  in  Phase  II  were: 

1)  VS2DTI  (USGS); 

2)  UNSATCHEM  (USSL); 

3)  HYDRUS-2D  (USSL-IGWMC), 

4)  SEVIEW(ESCI);  and, 

5)  CHEMFLOW2000  (OSU). 

Phase  II  evaluates  the  five  numerical  model  codes  identified  in  Phase  I  for  their  ability  to 
predict  the  shallow  movement  of  relatively  low  concentrations  of  salt  in  groundwater. 
Relatively  low  concentrations,  in  this  context,  is  interpreted  as  concentrations  below 
which  density-dependent  behaviour  can  be  ignored.  This  report  describes  the 
performance  of  the  five  codes  for  a  variety  of  generic  scenarios  with  varying  degrees  of 
complexity.  The  generic  scenarios  were  derived  from  a  review  of  typical  case  histories 
and  are  used  to  illustrate  and  compare  the  performance  of  each  of  the  five  selected 
computer  codes.  The  effects  of  density-dependent  flow  and  flow  through  discretely 
fractured  media  were  not  evaluated. 


A355-1 000002 
Page  1 


ENGINEERED  SOLUTIONS 


Evaluation  of  Computer  Models  -  Phase  H  Report 


March  2003 


The  principal  objective  of  Phase  II  of  the  project  is  to  document  the  strengths  and 
weaknesses  of  each  code  for  the  specific  task  of  simulating  the  transport  of  salt  (sodium 
chloride)  in  soil  and  groundwater  in  climatic  regimes  characteristic  of  the  Alberta  foothills 
and  prairies. 

2.0  SCOPE  OF  PHASE  II 

The  scope  of  the  Phase  II  project  was  to: 

1)  Evaluate  five  numerical  computer  codes  that  were  highly  ranked  in  the  pre-screening 
process  provided  in  Phase  I  (MDH,  2002). 

2)  Develop  generic  risked-based  scenarios  derived  from  typical  case  histories  to 
illustrate  movement  of  salt  from: 

a.  Leaks  in  shallow  buried  pipelines; 

b.  Localized  surface  sources  from  periodic  use  of  flare  pits;  and, 

c.  Long  term  surface  sources  from  highway  salt  storage  piles. 

3)  Provide  a  report  documenting: 

a.  The  information  requirements  to  model  each  scenario; 

b.  The  ability  of  each  code  to  simulate  the-  scenarios; 

c.  The  ease  of  use  of  each  code;  and 

d.  The  major  strengths  and  weaknesses  of  each  code. 

The  scope  of  Phase  II  did  not  examine  models  that  could  account  for  density  dependent 
flow  or  flow  through  fractured  media.  In  order  to  have  any  density  effect,  there  must  be 
a  density  contrast  between  the  contaminated  water  and  the  fresh  water.  That  means 
that  the  relative  density  of  the  contaminated  water  must  be  significantly  greater  or  less 
than  the  density  of  the  fresh  water.  If  there  is  no  significant  contrast  between  the 
contaminated  water  and  the  fresh  water,  then  the  problem  is  essentially  a  normal 
transport  problem  involving  advection  and  hydrodynamic  dispersion.  After  discussion 
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with  Alberta  Environment,  it  was  concluded  that  the  density  gradients  would  be  small  for 
the  scenarios  being  considered,  and  therefore,  density  dependent  flow  could  be 
neglected.  Consequently,  the  models  included  in  the  report  do  not  have  the  ability  to 
simulate  density  dependent  flow.  The  models  apply  to  circumstance  where  the 
concentration  of  salt  in  the  groundwater  does  not  create  large  density  gradients  that 
would  significantly  influence  the  long-term  movement  of  a  salt  plume.  If  the  site  being 
modelled  shows  signs  of  large  density  gradients,  more  complex  density-dependent 
codes  such  as  FEMWATER,  FEFLOW  or  SUTRA  should  be  used. 

The  effects  of  fracture  flow  were  not  included  in  the  study.  Clay-rich  deposits  cover 
large  areas  of  Canada  and  the  northern  United  States.  The  upper  layers  of  these 
deposits  are  typically  weathered  and  fractured.  Weathering  and  fracturing  can  increase 
the  otherwise  low  hydraulic  conductivity  of  the  deposits.  Grisak  and  Pickens  (1980) 
suggest  transport  of  dissolved  contaminants  in  a  fractured  clay  (or  porous  rock)  is 
expected  to  be  controlled  by  (1)  advection  through  the  fractures;  (2)  diffusion  into  the 
clay  matrix  between  fractures;  and  (3)  retardation  processes  in  both  the  fractures  and 
matrix.  There  are  only  a  few  field  experiments  that  have  estimated  the  hydraulic  fracture 
aperture  and  fracture  porosity  for  nonlithified  clay-rich  or  silt-rich  deposits.  Field  studies 
published  by  Keller  et.  al.  (1986)  and  McKay  et.  al.  (1993)  suggest  that  fractured  till 
could  have  a  representative  bulk  hydraulic  conductivity  of  up  to  3  orders  of  magnitude 
greater  than  the  mean  value  for  the  unfractured  matrix  determined  from  laboratory  tests. 
Often  the  fractures  in  the  till  have  been  incorporated  into  numerical  models  assuming  the 
fractures  are  frequent  and  pervasive  such  that  the  fractured  till  acts  as  an  "equivalent 
porous  medium"  or  epm,  and  can  be  accounted  for  by  increasing  the  hydraulic 
conductivity.  The  models  included  in  this  report  cannot  be  used  for  sites  that  contain 
media  where  the  individual  discrete  fractures  cannot  be  represented  as  part  of  the  soil 
matrix  in  terms  of  a  bulk  hydraulic  conductivity. 
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3.0  CODE  EVALUATION  FOR  RBCA 

The  primary  basis  of  the  evaluation  of  the  five  codes  is  their  suitability  for  screening  sites 
for  risk  based  corrective  action  (RBCA).  The  tiered  approach  of  RBCA  is  designed  to 
eliminate  low-risk  sites  from  further  consideration  with  minimum  analysis  and  expense 
while  retaining  high-risk  sites  for  progressively  more  detailed  review.  For  the  purposes 
of  this  report,  the  three  tiers  of  assessment  are  interpreted  as  follows: 

3.1  Tier  1  Assessment 

At  this  very  preliminary  level,  only  an  estimate  of  maximum  concentrations  and  total 
mass  at  the  point  of  exposure  (POE)  and  identification  of  receptors  are  required.  A 
simple  look  up  table  might  be  sufficient  to  decide  whether  the  maximum  concentration 
and/or  total  mass,  or  the  site  and/or  receptors,  justify  Tier  2  analysis.  Numerical 
modelling  is  not  usually  required  at  this  level. 

3. 2  Tier  2  Assessment 

At  this  level,  the  affected  porous  media  needs  to  be  delineated  together  with 
determination  of  the  site-specific  transport  processes  and  pathways  from  POE  to 
receptor.  Simple  1D  column  models  form  the  basis  for  the  decision  whether  or  not  to 
proceed  to  Tier  3  analysis. 

3.3  Tier  3  Assessment 

For  a  Tier  3  assessment,  the  affected  porous  media  needs  to  be  delineated  further 
together  with  more  thorough  determination  of  site-specific  transport  processes  and 
transport  pathways  from  POE  to  receptor.  More  complex  2D  or  3D  numerical  models 
can  be  applied  at  this  level  to  investigate  the  extent  of  the  problem  and  evaluate 
remediation  alternatives. 
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3.3.1    Information  Requirements 

The  same  general  information  is  required  to  complete  both  a  Tier  2  and  a  Tier  3 
analysis.  However,  more  detailed  site-specific  information  should  be  gathered  to 
complete  a  Tier  3  analysis. 

After  a  site  visit  by  a  professional  engineer  or  geoscientist,  the  following  information 
should  be  obtained  to  complete  Tier  2  analysis 

1)  A  professional  assessment  of  the  geology  to  a  depth  of  5  to  10  m  below  the  zone  of 
contamination; 

2)  A  professional  assessment  of  the  location  of  the  potentiometric  surface(s); 

3)  A  professional  assessment  the  temporal  and  spatial  extent  of  salt-sources;  and, 

4)  Site-specific  material  properties  based  on  literature  sources  and  engineering 
judgement. 

•  The  minimum  material  property  data  obtained  by  estimation  or  measurement 
should  include  bulk  density,  saturated  hydraulic  conductivity,  porosity,  water 
content  profile,  residual  water  content,  parameters  used  for  predicting  the  soil 
water  characteristic  curves  (SWCC),  and  dispersivities  and/or  dispersion 
coefficients. 

The  following  information  should  be  obtained  to  complete  Tier  3  analysis: 

1)  Borehole  information  indicating  geology  to  a  depth  of  5  to  10  m  below  the  zone  of 
contamination; 

2)  Borehole  information  adequate  to  construct  2D  site  cross-sections; 

3)  Location  of  the  potentiometric  surface(s)  and  groundwater  gradient(s); 

4)  Site  topography  and  aerial  photographs  to  locate  all  potential  receptors  near  the  site; 

•  If  aerial  photographs  or  topographic  maps  do  not  provide  appropriate  scale 
coverage  or  if  the  age  of  the  photographs  is  not  indicative  of  current  site 
conditions,  field  verification  would  be  necessary  to  ensure  receptors  are 
identified. 
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5)  Material  properties  based  on  site-specific  testing; 

•  The  minimum  material  property  data  should  include  bulk  density,  hydraulic 
conductivity,  porosity,  water  content  profile  and  residual  water  content. 

6)  Material  properties  based  on  literature  sources  and  engineering  judgement; 

•  Parameters  used  for  predicting  the  soil  water  characteristic  curves  (SWCC)  and 
dispersivities  and  diffusion  coefficients  could  be  measured,  but  are  generally 
taken  from  literature. 

7)  Site  history  sufficient  to  determine  the  temporal  and  spatial  extent  of  salt-sources; 
and, 

8)  Soil-salinity  data  indicating  contaminant  source  concentration  and/or  contaminant 
concentration  with  depth,  including  EM  surveys  and  groundwater  chemical  analyses. 

Gravels  and  sands  represent  high  hydraulic  conductivity  units  (aquifers)  and  tills,  clays, 
and  silts  comprise  the  low  hydraulic  conductivity  units  (aquitards).  Table  3.1  shows 
typical  ranges  of  hydraulic  conductivities  for  the  hydrostratigraphic  units  researched  from 
literature.  MDH  provides  these  values  only  as  a  rough  guide  and  strongly  advise  that 
modellers  should  use  their  own  judgement  in  the  application  of  these  values  to  the 
specific  site  being  modelled. 

Table  3.2  illustrates  values  of  porosity,  soil  compressibility  and  specific  storage 
researched  from  literature.  Specific  storage  is  the  volume  of  water  released  from  a  unit 
volume  of  confined  porous  medium  per  unit  decline  in  hydraulic  head  per  unit  thickness. 
Porosity  and  compressibility  of  both  water  and  porous  medium  are  related  to  specific 
storage  as  follows  (Freeze  and  Cherry,  1979): 


Ss  =  (a  +  np)Yv 


where:  Ss  =  specific  storage 

a  =  compressibility  of  the  soil  matrix 
P  =  compressibility  of  water 
n  =  porosity 

Yw  =  specific  weight  of  water 


[F-Y^] 
[] 


[FL-^] 
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Table  3.1  -  Hydraulic  conductivities  for  typical  hydrostratigraphic  units. 


Hydrostratigraphic  Units 

Hydraulic  Conductivity 

Upper  Limit 
(m/s) 

Lower  Limit 
(m/s) 

Surficia!  stratified  deposits 

sand 

A     A  n-7  (^) 

>1x10   ^  ' 

Sand 

5x1 0"^ 

n    A  n-7  (4)  (7) 

2x10   ^  ' 

Clay 

2x1 0'^ 

A  An-^^  (1)14) 
1x10    ^  ' 

Normally-consolidated  till 

Oxidized 

4x1 0'^ 

A    A  n-8  (^)  (7)  (8) 

1x10  ^''''''^°> 

Overconsolidated  till 

Oxidized 

4x1 0"^ 

A    An-8  (1)(7){8) 

1x10 

Unoxidized 

2x10-^ 

5^^Q-10  (1)(2)(7)(8) 

Intertill  glaciuofluvial  aquifers 

1x1 0'* 

3.5x10-^^^^ 

6x1 0'^ 

2x10"' 

Heavily  overconsolidated  clay  till 

1x10-^° 

^x^o■''  ^'^^'^ 

Buried  valley  aquifers 

6x10-^ 

2x10"' 

IxlO"* 

^x^Q■' 

Bedrock  sandstones 

1x10-^ 

1x10-^ 

Bedrock  shales 

1x10-^° 

Note:  The  references  apply  to  both  the  upper  and  lower  limits. 

1  Maathuis  and  van  der  Kamp  (1994a) 

2  Site  specific  data 

3  Keller  e^  a/.  (1988)  and  Keller  etal.  (1989) 

4  Domenico  and  Schwartz  (1998) 

5  Maathuis  and  Schreiner  (1982) 

6  Maathuis  ef  a/.  (1994) 

7  Ho  and  Barbour  (1987) 

8  Colder  (1997) 

9  Freeze  and  Cherry  (1979) 
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Table  3.2  -  Storage  properties  for  typical  hydrostratigraphic  units. 


Hydrostratigraphic  Unit 

Compressibility 
(m'/N) 

Porosity 

(%) 

Specific  Storage 

Upper 
Range 

Lower 
Range 

Upper 
Range 

Lower 
Range 

Upper  Range 

Lower 
Range 

Surficial  stratified 
deposits 

sand 

1.5x10'^ 

1.5x10"^ 

50 

25 

1.9x10"^ 

6.9x10^'^' 

silt 

1.2x10"^ 

1.2x10"^ 

50 

35 

1 .9x1 0"^ 

7.4x10  *  ' 

clay 

1.5x10"^ 

1.5x10"^ 

50 

40 

2.7x10'^ 

9.5x10"^ 

Normally  consolidated  till 

6.9x10"^ 

1.3x10"^ 

40 

20 

1.3x10"^ 

2.8x1 0"^<^' 

Overconsolidated  till 

6.9x10"^ 

1.3x10"® 

45 

20 

1.3x10"^ 

2.8x10"^^^' 

Intertill  aquifer 

1.5x10"^ 

1.5x10"^ 

30 

25 

3.8x10"^ 

1.7x10"' 

Heavily  overconsolidated  till 

1.2x10-^ 

1.2x10"® 

40 

30 

2.7x10"^ 

5.9x10"^ 

Buried  valley  aquifer 

1.5x10"^ 

5.5x10"^ 

35 

20 

1.0x10"^ 

3.8x10"^ 

Bedrock  sandstone 

1.0x10-^ 

I.OxlO"""" 

20 

5 

1.0x10"^ 

1.0x10"^ 

Bedrock  shale 

1.0x10"^ 

I.OxlO"''" 

15 

5 

1.0x10"^ 

1.0x10"^ 

Notes:    Compressibility  of  water  at  25°C  is  4.8  x  1 0""""  m^/N 

Soil  compressibility  values  estimated  from  Dominico  and  Schwartz  (1990)  and  Canmet  (1986) 

1  Maathuis  and  van  der  Kamp  (1994b) 

2  Therrien  and  Sudicky  (1996) 

3  Freeze  and  Cherry  (1979) 


For  unconsolidated  sediments,  the  compressibility  of  the  soil  is  much  greater  than  that  of 
water  so  the  specific  storage  is  controlled  by  the  compression  of  the  pore  skeleton.  For 
rigid  bedrock  aquifers,  the  opposite  is  true  and  the  specific  storage  is  controlled  by  the 
compressibility  of  the  pore  fluid.  A  reduction  in  head  or  pore-pressure  results  in  a 
release  of  fluid  from  storage  and  a  reduction  in  storage  volume  as  the  skeleton 
collapses.  An  increase  in  head  or  pore-pressure  corresponds  to  an  expansion  of  the 
storage  volume  as  the  skeleton  is  forced  to  dilate. 

Porosities  tend  to  lie  in  a  relatively  narrow  range  from  5  to  50%  (Fetter,  1994; 
Freeze  &  Cherry,  1979).  Porosity  values  can  be  estimated  based  on  typical  ranges  for 
different  lithologies  as  cited  in  the  literature.  Porosity  is  a  key  parameter  in  contaminant 
migration  and  fractured  material  can  exhibit  what  is  called  "dual  porosity".  Clay  tills  can 
have  matrix  porosities  up  to  50%  with  a  superimposed  fracture-porosity  of  less  than  5%. 
Using  the  matrix  porosity  rather  than  the  fracture  porosity  in  contaminant  transport 

^MDH 
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calculations  leads  to  an  initial  underestimate  of  the  rate  of  movement  since  transport 
takes  place  rapidly  through  the  fracture  network  and  displaces  the  matrix  pore  fluids 
more  slowly.  However,  none  of  the  evaluated  codes  have  the  capability  to  model  "dual 
porosity"  materials. 

The  hydraulic  conductivities  of  soils  in  the  unsaturated  zone  are  lower  than  the  same  soil 
under  saturated  conditions  (below  the  water  table).  The  reduced  moisture  contents  in 
the  unsaturated  zone  result  in  less  of  the  pore-network  being  available  for  flow.  Suction 
pressures  (negative  pressure  heads)  are  preset  in  the  unsaturated  zone,  with  the  water 
table  (matric  suction  =  0)  marking  the  intert'ace  between  saturated  and  unsaturated 
conditions.  Hydraulic  conductivity  is  highly  dependent  on  suction  pressure  when  a  soil  is 
unsaturated  and  sands  and  gravels  in  the  unsaturated  state  can  be  barriers  to  flow. 
Values  for  typical  coarse  and  fined  grained  materials  are  presented  in  Table  3.3. 


Table  3.3  -  Van  Genuchten  Parameters  for  prediction  of  SWCC. 


Material 

Alpha  (m"') 

Beta 

Fine  Grained  (clay  till) 

1.5 

1.1 

Medium  Grained  (silt,  silty  clay) 

5 

1.8 

Coarse  Grained  (sand) 

12 

2.7 

1  Freeze  and  Cherry  (1 979) 

2  Domenico  and  Schwartz  (1998) 


Typical  pressure  head  versus  relative  hydraulic  conductivity  (K(0)/Ksat),  degree  of 
saturation  (Se),  and  moisture  content  (0)  curves  for  coarse  and  fine  grained  materials  are 
provided  in  Figure  3.1  (a)  and  (b),  respectively. 
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> 

■g  1.0E-03 


50  75 
Saturation  (%) 


(a)  Coarse  grained  material 


Saturation  (%) 


94  96 
Saturation  (%) 


(b)  Fine  grained  material 
Figure  3.1  -  Sample  soil  water  characteristic  curves. 


25  50  75 

Saturation  (%) 
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Sand  is  easily  drained  showing  a  rapid  change  in  saturation  for  a  small  suction  head 
increment  and  a  corresponding  rapid  drop  in  the  relative  hydraulic  conductivity.  The 
fine-grained  material  remains  almost  fully  saturated  even  at  large  suction  heads  and  only 
a  modest  change  in  relative  hydraulic  conductivity  is  associated  with  drainage.  The 
curves  were  developed  using  Van  Genuchten  parameters  based  on  fitting  the  results  of 
laboratory  tests  on  similar  materials.  The  program  SoilVision  (SoilVision  Systems  Ltd., 
1997)  and  its  associated  database  were  used  to  assist  in  this  process. 

Van  Genuchten  (1980)  developed  an  empirical  formula  relating  hydraulic  conductivity 
and  moisture  content  with  negative  pressure  head: 


where:  6  =  volumetric  moisture  content 
6r=  residual  moisture  content 
0s  =  saturated  moisture  content 
a  =  van  Genuchten  curve-fitting  parameter  [L'^] 
\|/  =  matric  potential  or  suction  pressure  [L] 
p  =  van  Genuchten  curve-fitting  parameter 
m  =  1-1/p 

Van  Genuchten's  function  was  developed  from  Mualem's  (1978)  relationship,  used  to 
predict  the  hydraulic  conductivity  from  moisture  retention  data: 


where:  Se  =  (0  -  0r  )/(0s  -  0r )  is  the  degree  of  saturation 

L  =  is  a  pore-conductivity  parameter  or  tortuosisy  factor 

Using  Mualem's  (1978)  model  Van  Genuchten  (1980)  derived  a  closed-form  analytical 
solution  to  predict  the  relative  hydraulic  conductivity  (Kr)  at  a  given  volumetric  water- 
content: 


0  =  0r  +  ( 0.  -  0r ) 
[1+laiKp] 


m 


Se 


K(Se)  =  KsatSe'  [f(Se)  /  f  (1 )]'  where  f  (Se)  =  J  1  /  v|/(x)  dx 


0 


Kr  =  Se'{1  -[1  -Se'^"]  ""f 
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Where:  Se  can  also  be  written  as  [1  +  |  a\]/|  ^ 

The  confined  storage  characteristics  of  the  hydrostratigraphic  units  can  be  determined 
by  multiplying  specific  storage  by  the  thickness.  On  the  water  capacity  against  depth 
curves  the  specific  storage  corresponds  to  the  gradient  at  the  water  table. 

4.0  GENERIC  SCENARIOS 

The  first  set  of  generic  scenarios  comprised  simple  1D  columns  with  a  (non-reactive  or 
conservative)  tracer  contaminant.  Analysis  of  the  1D  columns  was  completed  to  verify 
that  CHEMFLO,  HYDRUS,  SEVIEW,  UNSATCHEM,  and  VS2DTI  performed 
consistently  using  the  same  input  data.  The  initial  comparison  was  completed  using 
generic  data  for  both  homogeneous  (uniform)  and  heterogeneous  (layered)  systems 
including  both  fine  and  coarse-textured  soils. 

Additional  2D  radial  (axisymmetric)  analysis  was  then  completed  using  HYDRUS  and 
VS2DTI  on  generic  scenarios  involving  more  complex  sources  and  boundary  conditions. 

Finally  scenarios  involving  cation-exchange  reactions  using  UNSATCHEM  and  VS2DTI 
were  completed  to  evaluate  the  ability  of  these  codes  to  simulate  simple  reactive 
transport  scenarios  involving  cation  exchange. 

4.1  1D  Column  Scenarios 

This  phase  of  modelling  was  carried  out  to  test  the  ability  of  the  codes  to  produce  the 
same  results  using  equivalent  or  near-identical  input.  Initially  generic  simulations  were 
carried  out  on  simple  1D  columns  with  a  homogeneous  isotropic  coarse  or  fine-grained 
soil  to  verify  that  all  five  codes  produce  consistent  results.  For  these  analyses  it  was 
anticipated  that  the  results  would  be  the  same  because  the  geology  and  mesh- 
discretization  are  straightfonA/ard  and  the  boundary  conditions  are  simple  and 
unambiguous. 

The  next  step  involved  increasing  the  complexity  of  the  1D  columns  to  include  three 
layers  of  varying  lithology  with  both  saturated  and  unsaturated  flow  zones.  Such 
scenarios  represent  typical  Tier  2  screening  tasks  where  the  local  stratigraphy  can  be 
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deduced  and  generic  material-properties  can  be  assigned  by  an  experienced 
practitioner.  From  this  phase  of  the  analysis,  an  initial  comparison  of  the  codes  was 
made  in  terms  of  ease-of-use,  flexibility,  convergence  characteristics,  and  mass  balance 
performance. 

All  columns  were  30  m  high  (deep)  with  a  nominal  cross  section  of  1  m^.  The  nominal 
dimensions  of  the  columns,  material  properties,  solute  concentrations  and  application 
times  were  kept  constant  for  each  simulation.  The  surface-flux  boundary  condition  was 
changed  depending  on  whether  a  coarse  or  fine-grained  soil  appeared  at  the  surface  to 
simulate  different  rates  of  recharge.  In  every  case  the  water  table  was  assumed  to  be  at 
a  constant  depth,  arbitrarily  chosen  to  be  3  m  below  ground  surface  (bgs). 

All  heterogeneous  columns  were  initially  run  as  "flow  only"  problems  to  obtain  an 
approximate  steady-state  initial  hydraulic  condition  for  the  transport  simulation.  After  a 
steady-state  solution  was  obtained,  a  source  with  a  concentration  corresponding  to  brine 
with  a  density  similar  to  seawater  (30  kg/m^,  30  g/L  or  30,000  mg/L)  was  applied  to  the 
surface  of  the  column  for  20  years.  After  20  years,  the  source  was  removed  and  the 
model  was  run  for  an  additional  20  years  to  simulate  movement  of  the  distributed 
source. 

The  source  was  applied  as  a  constant  concentration  and  constant  flux  boundary.  This 
boundary  condition  was  chosen  because  the  salt  was  assumed  to  enter  the  column  at  a 
rate  determined  by  the  ability  of  infiltration  to  enter  the  soil.  Neither  constant 
concentration  nor  constant  head  boundaries  were  used  because  these  boundaries 
generally  apply  to  a  ponded  source  where  the  inflow  rate  is  determined  by  the  saturated 
hydraulic  conductivity  of  the  soil.  This  was  not  likely  to  be  the  case  for  a  shallow  pipe 
break,  intermittent  flare  pit  or  "dry"  salt  pile.  Although  the  water  may  be  ponded  in  a  flare 
pit,  it  was  assumed  that  the  ponding  would  not  occur  for  long  enough  to  maintain  the 
water  table  at  the  surface  for  an  extended  period. 

For  1D  columns,  it  is  necessary  to  specify  a  second  boundary  condition  at  the  base  of 
the  column.  Such  boundaries  may  be  specified  as  free  draining  or  as  an  advective  flow 
with  a  user-specified  concentration.  A  third  choice  is  a  constant  head  condition  where 
the  flow  and  mass  flux  are  computed  by  the  code.  For  the  comparative  simulations,  a 
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constant  pressure  head  boundary  was  applied  to  the  base  of  the  1D  columns 
corresponding  to  the  height  of  the  water  table. 

The  material  properties  for  the  generic  coarse  and  fine-grained  soils  are  shown  in 
Table  4.1  and  Table  4.2,  respectively. 


Table  4.1  -  Material  properties  for  a  homogeneous  isotropic  coarse-grained  soil 

(sand). 


Input  Parameter 

Value 

Units 

Dimension 

Hydraulic  Conductivity 

6.9 

m/d 

L/T 

Saturated  Water  Content  (Porosity) 

0.25 

fraction 

dimensionless 

Residual  Water  Content 

0.05 

fraction 

dimensionless 

Van  Genuchten  parameter  alpha 

12 

1/m 

1/L 

Van  Genuchten  parameter  beta 

2.7 

exponent 

dimensionless 

The  material  properties  in  Table  4.1  were  chosen  to  be  characteristic  of  a  clean  medium- 
grained  sand.  Such  sand  horizons  commonly  occur  in  glaciofluvial  channels. 


Table  4.2  -  Material  properties  for  a  homogeneous  isotropic  fine-grained  soil  (clay 

till). 


Input  Parameter 

Value 

Units 

Dimension 

Hydraulic  Conductivity 

8.6  X  10"^ 

m/d 

L/T 

Saturated  Water  Content  (Porosity) 

0.36 

fraction 

dimensionless 

Residual  Water  Content 

0.30 

fraction 

dimensionless 

Van  Genuchten  parameter  alpha 

1.5 

1/m 

1/L 

Van  Genuchten  parameter  beta 

1.1 

exponent 

dimensionless 

The  materials  in  Table  4.2  were  chosen  to  be  characteristic  of  a  poorly-sorted 
homogeneous  till  with  a  silty-clay  matrix.  It  is  assumed  that  the  till  is  a  single-porosity 
porous  medium.  A  dual  porosity  medium  with  a  component  of  fracture  flow  is  not 
considered.  For  such  a  material,  a  considerably  lower  porosity  might  be  used  to 
characterize  the  solute  interaction  with  the  fracture-porosity  component  alone. 
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The  solute  transport  parameters  used  in  the  generic  models  are  listed  in  Table  4.3. 
Diffusion  coefficients  for  porous  media  were  estimated  by  multiplying  the  self-diffusion 
coefficient  by  porosity.  A  single  (low)  value  was  applied  for  longitudinal  dispersivity. 
Dispersivity  is  generally  regarded  as  a  scale-dependent  property  and  the  value  used  is 
appropriate  on  the  tens-of-metres  scale.  For  solute  transport  on  a  scale  of  hundreds  of 
metres,  a  second  value  is  included  in  the  table. 


Table  4.3  -  Solute  transport  parameters. 


Input  Parameter 

Value 

Units 

Dimension 

Diffusion  Coefficient  (CI  ) 

4.4  X  10-^ 

m^/d 

L^/T 

Diffusion  Coefficient  (Na^) 

2.9  X  10-^ 

m'/d 

L^/T 

Diffusion  Coefficient  (Ca^"") 

1.7  X  10-^ 

m^/d 

Longitudinal  Dispersivity  (10  m  scale) 

0.4 

m 

L 

Longitudinal  Dispersivity  (100  m  scale) 

4.0 

m 

L 

Dispersivity  is  the  property  of  a  porous  medium  that  characterizes  mechanical  mixing  in 
a  flow  field.  Longitudinal  dispersivity  (aL)  quantifies  mixing  in  the  direction  of  flow  and 
transverse  dispersivity  (aj)  characterizes  mixing  normal  to  the  mean  flow  direction.  In 
common  with  hydraulic  conductivity,  dispersivity  varies  with  the  scale  of  measurement. 
A  reasonable  rule  of  thumb  is  that  macroscopic  dispersion  (field  scale)  is  approximately 
2  orders  of  magnitude  greater  than  microscopic  dispersion  (lab  scale).  At  the  lab  column 
scale,  longitudinal  dispersivity  is  typically  between  0.01  and  1.0  cm.  In  field  studies, 
values  from  0.1  to  2.0  m  are  observed  over  short  distances  in  well-controlled 
experiments.  Longitudinal  dispersivity  values  >10  m  are  reported  over  longer  distances, 
but  data  is  very  sparse  and  based  on  the  back-analysis  of  documented  plumes. 

The  famous  Borden  experiment  (Sudickey,  1986;  Mackay  et  ai,  1986)  conducted  field 
scale  dispersivity  measurements.  The  study  found  that  medium-grained,  fine-grained, 
and  silty  fine-grained  sand,  gave  a  field  dispersivity  value  of  0.45  m  over  a  travel 
distance  58  m.  Transverse  dispersivity  is  generally  at  least  an  order  of  magnitude 
smaller  than  longitudinal  dispersivity.  Table  4.4  summarizes  the  information  used  to 
estimate  dispersivities. 
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Table  4.4  -  Solute  transport  parameters. 


Dispersivity  Scale 

Value  Range  (m) 

Laboratory  homogeneous  sand  column  experiments 

10  "^  to  10"' 

Naural-gradient  tracer  field  experiments 

10"'  to  10"° 

Empirical  fitting  or  matching  of  field  plumes 

10°  to  10' 

It  is  well  established  that  the  apparent  dispersivity,  used  to  compute  the  dispersion 
coefficient  in  the  advection-dispersion  equation,  exhibits  scale  effects.  The  apparent 
dispersivity  grows  with  the  plume  size  until  the  size  of  the  plume  reaches  10-15  times 
the  size  of  the  largest  scale  of  heterogeneity  (Howington  et  al.,  1997).  For  typical  sand 
aquifers,  interpreting  the  largest  heterogeneity  as  lenticular  bodies  with  a  length  scale  of 
1-5  m,  an  appropriate  choice  for  dispersivity  would  seem  to  be  about  10  to  50  m. 

The  information  common  to  all  the  codes  is  provided  in  Tables  4.1  to  4.3.  However, 
each  code  requires  the  input  data  in  a  slightly  different  form.  Therefore,  a  detailed 
description  of  the  information  specific  to  each  modelling  code  is  provided  in  Appendix  A. 

4.1.1    Homogeneous  Isotropic  Coarse-Grained  Soil  Column  -  Surface  Source 

The  boundary  conditions  utilized  in  the  modelling  of  the  uniform  isotropic  coarse-grained 
soil  column  are  provided  in  Table  4.5.  An  exit-solution  concentration  of  zero  was  applied 
to  the  base  of  the  column.  It  was  assumed  that  the  base  of  the  column  was  below  the 
zone  of  contamination,  and  no  contaminant  would  exit  in  the  outflow  during  the 
simulation  period. 


Table  4.5  -  Boundary  conditions  for  a  homogeneous  isotropic  coarse-grained  soil 

column  (surface  source). 


Boundary 

Boundary  Type 

Value 

Units 

Surface  Boundary  for  Flow 

Constant  Flux 

5.4  x  10"^ 

m/d 

Base  Boundary  for  Flow 

Constant  Head 

27 

m 

Surface  Boundary  for  Solute  (0-20  years) 

Solute  Concentration 

30 

kg/m^ 

Surface  Boundary  for  Solute  (20-40  years) 

Solute  Concentration 

0 

kg/m^ 
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The  surface  flux  is  characteristic  of  infiltration  rates  into  sandy  surfaces  under  climatic 
conditions  expected  in  the  western  Canadian  prairies.  Figure  4.1  illustrates  the 
concentration  distribution  of  the  contaminant  plume  after  20  years  of  source  application 
for  CHEMFLO,  HYDRUS,  UNSATCHEM,  and  VS2DI. 


^  10 
E 


£  15 

Q. 

0) 

°  20 

25 
30 


CIC. 


■CHEMFLO 


VS2DTI 


HYDRUS 


UNSATCHEM 


Figure  4.1  -  Concentration  versus  depth  profiles  after  20  years  of  source 
application  (homogeneous  isotropic  coarse-grained  soil  column). 


Figure  4.1  shows  that  the  modelling  results  for  CHEMFLO,  HYDRUS,  UNSATCHEM  and 
VS2DTI  correspond  very  closely.  For  the  coarse-grained  soil  column,  the  plume  has 
penetrated  to  a  depth  of  approximately  6.8  metres  (from  6.7  to  6.8  m)  after  20  years  of 
continuous  source  applied  to  the  surface.  The  peak  concentration  is  approaching  the 
source  concentration  of  30  kg/m^.  Plume  penetration  depth  is  arbitrarily  defined  as  the 
depth  where  concentration  exceeds  1  kg/m^  or  C/Co  >  0.033.  Detailed  input  and  output 
data  for  the  simulations  can  be  found  in  Appendix  A  and  Appendix  B,  respectively. 

Results  are  not  shown  for  SEVIEW  (SESOIL  combined  with  ATI 23).  SESOIL  computes 
a  water  balance  from  seasonal  climatic  data  and  applies  the  groundwater  recharge  to 
the  surface  of  the  ATI 23  flow  model.  SESOIL  is  designed  as  a  screening  program  for 
low  concentration  organic  contaminants  in  relatively  permeable  surface  soils.  The 
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application  of  SESOIL  and  AT123  are  limited  to  homogeneous  materials  and  can  only 
provide  a  Tier  1  level  of  screening.  SESOIL  was  unsuccessful  in  predicting  positive 
recharge  for  raw  climatic  data  from  the  western  Canadian  prairies.  After  considerable 
data  manipulation,  it  was  possible  to  generate  positive  groundwater  recharge  but  it 
proved  extremely  difficult  to  apply  SESOIL  to  commonly  occurring  low-permeability  soils 
(e.g.  tills).  For  these  reasons,  comparative  evaluation  of  SEVIEW  and  the  other  codes 
was  not  feasible. 

4.1.2    Homogeneous  Isotropic  Fine  Grained  Soil  Column  -  Surface  Source 

The  boundary  conditions  utilized  in  modelling  the  homogeneous  fine-grained  soil  column 
are  provided  in  Table  4.6.  An  exit-solution  concentration  of  zero  was  applied  to  the  base 
of  the  column,  because  it  was  assumed  that  the  base  of  the  column  was  below  the  zone 
of  contamination  and  no  contaminant  would  exit  in  the  outflow  during  the  simulation 
period. 


Table  4.6  -  Boundary  conditions  for  a  homogeneous  isotropic  fine  grained  soil 

column  (surface  source). 


Boundary 

Boundary  Type 

Value 

Units 

Surface  Boundary  for  Flow 

Constant  Flux 

1.0x10"^ 

m/d 

Base  Boundary  for  Flow 

Constant  Head 

27 

m 

Surface  Boundary  for  Solute  (0-20  years) 

Solute  Concentration 

30 

kg/m^ 

Surface  Boundary  for  Solute  (20-40  years) 

Solute  Concentration 

0 

kg/m^ 

The  surface  flux  value  is  considered  characteristic  of  infiltration  rates  into  till  surfaces 
under  climatic  conditions  expected  in  the  western  Canadian  prairies.  Figure  4.2 
illustrates  the  concentration  of  the  contaminant  plume  after  20  years  of  source 
application. 

Figure  4.2  shows  that  the  modelling  results  for  CHEMFLO,  HYDRUS,  UNSATCHEM  and 
VS2DTI  correspond  very  closely.  For  the  fine-grained  soil  column,  the  plume  has 
penetrated  to  a  depth  of  approximately  2.0  metres  (1.8  to  2.1  m)  after  20  years  of 
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continuous  source  applied  to  the  surface.  The  C/C°  is  approximately  0.36.  Detailed 
input  and  output  data  is  provided  in  Appendix  A  and  Appendix  B,  respectively. 
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Figure  4.2  -  Concentration  versus  depth  profiles  after  20  years  of  source 
application  (homogeneous  isotropic  fine-grained  soil  column). 


4.1.3    Heterogeneous  (Layered  Fine-Coarse-Fine)  Isotropic  Soil  Column  -  Surface 
Source 

The  boundary  conditions  utilized  in  the  modelling  of  the  isotropic,  layered  fine-coarse- 
fine  soil  column  are  provided  in  Table  4.7.  An  exit-solution  concentration  of  zero  was 
applied  to  the  base  of  the  column,  because  it  was  assumed  that  the  base  of  the  column 
was  below  the  zone  of  contamination,  and  no  contaminant  would  exit  in  the  outflow 
during  the  simulation  period. 
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Table  4.7  -  Boundary  conditions  for  an  isotropic,  layered  fine-coarse-fine  soil 

column  (surface  source). 


Boundary 

Boundary  Type 

Value 

Units 

Surface  Boundary  for  Flow 

Constant  Flux 

1.0  X  10"^ 

m/d 

Base  Boundary  for  Flow 

Constant  Head 

27 

m 

Surface  Boundary  for  Solute  (0-20  years) 

Solute  Concentration 

30 

kg/m^ 

Surface  Boundary  for  Solute  (20-40  years) 

Solute  Concentration 

0 

kg/m^ 

The  surface  flux  utilized  in  the  models  is  characteristic  of  infiltration  rates  into  till  surfaces 
under  clinnatic  conditions  expected  in  the  western  Canadian  prairies.  Figure  4.3 
illustrates  the  concentration  of  the  contaminant  plume  after  20  years  of  source 
application. 


Figure  4.3  -  Concentration  versus  depth  profiles  after  20  years  of  source 
application  (layered  fine-coarse-fine  soil  column). 

Figure  4.3  shows  that  the  modelling  results  for  CHEMFLO,  HYDRUS,  UNSATCHEM  and 
VS2DTI  correspond  very  closely.  For  the  fine-coarse-fine  soil  column,  the  plume  has 
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penetrated  to  a  depth  of  approximately  2.0  m  (from  1.8  to  2.1  m)  after  20  years  of 
continuous  source  applied  to  the  surface.  The  C/C°  is  about  0.36.  The  result  is  very 
similar  to  the  uniform  fine-grained  soil  because  the  plume  has  not  reached  the  first  fine- 
coarse  interface  (at  a  depth  of  4  m)  after  the  20  years  of  simulation.  Detailed  input  and 
output  data  is  provided  in  Appendix  A  and  Appendix  B,  respectively. 

4.1.4    Heterogeneous  (Layered  Coarse-Fine-Coarse)  Isotropic  Soil  Column  - 
Surface  Source 

The  boundary  conditions  utilized  in  the  modelling  of  the  isotropic,  layered  coarse-fine- 
coarse  soil  column  with  a  surface  source  are  provided  in  Table  4.8.  An  exit-solution 
concentration  of  zero  was  applied  to  the  base  of  the  column,  because  it  was  assumed 
that  the  base  of  the  column  was  below  the  zone  of  contamination,  and  no  contaminant 
would  exit  in  the  outflow  during  the  simulation  period. 


Table  4.8  -  Boundary  conditions  for  an  isotropic,  layered  coarse-fine-coarse  soil 

column  (surface  source). 


Boundary 

Boundary  Type 

Value 

Units 

Surface  Boundary  for  Flow 

Constant  Flux 

5.4x10'^ 

m/d 

Base  Boundary  for  Flow 

Constant  Head 

27 

m 

Surface  Boundary  for  Solute  (0-20  years) 

Solute  Concentration 

30 

kg/m^ 

Surface  Boundary  for  Solute  (20-40  years) 

Solute  Concentration 

0 

kg/m^ 

The  surface  flux  applied  to  the  soil  column  is  characteristic  of  infiltration  rates  into  sandy 
surfaces  under  climatic  conditions  expected  in  the  western  Canadian  prairies. 
Figure  4.4  illustrates  the  concentration  of  the  contaminant  plume  after  20  years  of  source 
application. 

Figure  4.4  shows  that  the  modelling  results  for  CHEMFLO,  HYDRUS,  UNSATCHEM  and 
VS2DTI  correspond  very  closely.  For  the  coarse-fine-coarse  soil  column,  the  plume  has 
penetrated  to  a  depth  of  about  6.2  m  (from  6.1  to  6.2  m)  after  20  years  of  continuous 
source  applied  to  the  surface.    The  peak  concentration  is  close  to  the  source 
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concentration  at  about  30  kg/m^.  The  plume  penetration  depth  is  approximately  0.5  m 
higher  than  the  uniform  coarse  soil  case  because  the  plume  is  slowed  at  the  first  coarse- 
fine  interface  at  a  depth  of  4  m.  Detailed  input  and  output  data  is  provided  in  Appendix  A 
and  Appendix  B,  respectively. 
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Figure  4.4  -  Concentration  versus  depth  profiles  after  20  years  of  source 
application  (layered  coarse-fine-coarse  soil  column). 

Figure  4.5  shows  the  progress  of  the  contaminant  plume  over  a  period  of  40  years  (20 
years  of  source  application).  Profiles  after  5,  10,  15,  20,  25,  30,  35  and  40  years  are 
provided  in  Figure  4.5.  The  results  are  plotted  for  UNSATCHEM  but  are  very  similar  for 
all  codes.  Over  the  period  of  source  application,  the  leading  edge  of  the  plume  sinks 
from  surface  to  6.2  m.  After  the  source  is  removed,  the  profiles  show  the  broadening 
and  attenuation  of  the  plume  peak  over  time  from  an  initial  value  of  30  kg/m^  (or  C/Co  = 
1.00  at  surface)  to  13  kg/m^  (or  C/Co  =  0.43  at  4.3  m  depth  after  40  years).  The  leading 
edge  of  the  plume  has  reached  a  depth  of  8.6  m  after  40  years  of  simulation. 
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Figure  4.5  -  UNSATCHEM  concentration  versus  deptli  profiles  (layered  coarse- 
fine-coarse  soil  column  with  20  year  source  application). 


4.1.5    Isotropic  Layered  Coarse-Fine-Coarse  Soil  Column  -  Distributed  Source 

The  boundary  conditions  utilized  in  the  modelling  of  the  isotropic,  layered  coarse-fine- 
coarse  soil  column  are  provided  in  Table  4.9.  An  exit-solution  concentration  of  zero  was 
applied  to  the  base  of  the  column  because  it  was  assumed  that  the  base  of  the  column 
was  below  the  zone  of  contamination,  and  no  contaminant  would  exit  in  the  outflow 
during  the  simulation  period.  The  distributed  source  varied  linearly  from  20  kg/m^  at 
surface  to  30  kg/m^  at  a  depth  of  3  m.  The  source  then  reduced  linearly  to  zero  at  a 
depth  of  10  m. 
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Table  4.9  -  Boundary  conditions  for  isotropic  layered  course-fine-course  soil 
column  (distributed  source). 


Dountiary 

Duuiiudry  1  ypc 

VdlUc 

Unite 

Surface  Boundary  for  Flow 

Constant  Flux 

5.4  X  10'^ 

m/d 

Base  Boundary  for  Flow 

Constant  Head 

27 

m 

Surface  Boundary  for  Solute  (0  years)* 

Solute  Concentration 

20 

kg/m^ 

3  m  Depth  Boundary  (0  years)* 

Solute  Concentration 

30 

kg/m^ 

10  m  Depth  Boundary  (0  years)* 

Solute  Concentration 

0 

kg/m^ 

Surface  Boundary  for  Solute  (20-40  years) 

Solute  Concentration 

0 

kg/m^ 

*  Note:  Initial  concentrations  varied  linearly  between  depths. 


The  surface  flux  is  characteristic  of  infiltration  rates  into  sandy  surfaces  under  climatic 
conditions  expected  in  the  western  Canadian  prairies.  Figure  4.6  illustrates  the 
concentration  of  the  contaminant  plume  after  20  years. 
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After  20  years  of  simulation,  the  source  distributions  predicted  by  CHEMFLO,  HYDRUS, 
UNSATCHEM  and  VS2DTI  were  all  similar,  with  a  peak  relative  concentration  (C/Co)  of 
0.77  (from  0.73  to  0.83)  at  a  depth  of  5.2  metres  (5.1  to  6.5  m).  CHEMFLO  appears  to 
predict  slightly  faster  migration,  with  the  peak  concentration  deeper  than  predicted  using 
the  other  codes  (Figure  4.6).  This  is  likely  due  to  small  differences  in  the  location  of  the 
distributed  source  relative  to  the  model  mesh  or  grid  in  the  different  simulations. 
Detailed  input  and  output  data  is  provided  in  Appendix  A  and  Appendix  B,  respectively. 

4.2  2D  Axisymmetric  (Radial)  Scenarios 

After  the  initial  comparison  and  verification  of  consistency  was  completed  on  the  1D  soil 
columns,  a  comparison  was  made  of  the  ability  of  the  codes  to  provide  a  Tier  3  site 
assessment.  The  pipeline,  flare  pit,  and  salt  pile  scenarios  were  simulated  using  2D 
axisymmetric  models  with  three  material  layers.  Of  the  five  codes,  only  HYDRUS-2D 
and  VS2DTI  were  capable  of  this  level  of  analysis.  Some  variation  in  results  between 
codes  was  anticipated  for  these  scenarios  because  each  program  has  different  input 
parameters  and  specification  details  for  mesh  construction,  boundary  conditions  and 
material  properties. 

4.2.1    Pipeline  Breal^  Scenario 

The  pipeline  scenario  was  modeled  for  an  isotropic,  heterogeneous  (coarse-fine-coarse 
layered)  soil  section.  The  pipe  (with  a  nominal  diameter  of  500  mm)  was  placed  with  the 
centreline  at  a  depth  of  1.25  m  below  ground  surface.  The  pipeline  scenario  was 
modelled  assuming  the  pipe  break  to  be  a  finite  length  (500  mm),  with  a  constant 
leakage  rate  and  duration.  Because  of  the  limitations  of  an  axisymmetric  model,  the 
pipe  break  was  limited  to  being  modeled  as  a  point  source.  The  source  appears  as  a 
"cylinder"  in  the  model,  with  a  height  and  diameter  of  500  mm  as  shown  in  Figures  4.7 
and  4.8. 
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Figure  4.7  -  2D  schematic  of  pipeline  breal<  (red  blocl^  indicates  breal<  location). 


Figure  4.8  -  3D  radial  schematic  of  pipeline  break  (red  block  indicates  break 

location). 

The  boundary  conditions  utilized  in  modelling  the  pipeline  break  in  the  layered  coarse- 
fine-coarse  soil  are  provided  in  Table  4.10. 


Table  4.10  -  Boundary  conditions  for  pipeline  break  scenario. 


Boundary 

Boundary  Type 

Value 

Units 

Surface  Boundary  for  Flow 

Constant  Flux 

5.4  X  10-^ 

m/d 

Radial  Boundary  for  Flow  -  Upper  Aquifer 

Constant  Total  Head 

-3 

m 

Radial  Boundary  for  Flow  -  Lower  Aquifer 

Constant  Total  Head 

-4 

m 

Flow  from  Pipe  into  System 

Constant  Flow 

0.02 

m^/d 

Pipe  Boundary  for  Solute  (0  - 1  years) 

Solution  Concentration 

30 

kg/m^ 

Pipe  Boundary  for  Solute  (1-11  years) 

Solution  Concentration 

0 

kg/m^ 
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Figure  4.9  illustrates  that  after  1  year  of  leakage,  HYDRUS  and  VS2DTI  produce  similar 
predictions  for  the  profile  and  relative  concentration  of  the  plume.  The  peak 
concentration  beneath  the  pipeline  break  occurs  at  a  depth  of  approximately  2.5  m. 
Both  codes  predict  salt  concentrations  of  about  5  to10%  of  the  source  at  the  ground 
surface. 
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Figure  4.9  -  Concentration  versus  depth  for  pipeline  break  after  1  year. 


Figure  4.10  (HYDRUS)  and  Figure  4.11  (VS2DTI)  show  that  the  plume  extends  to  a 
radial  distance  of  approximately  3.0  m  after  1  year  of  leakage.  After  this  time,  the  source 
concentration  was  turned  off  and  the  simulation  was  run  for  another  10  years.  Detailed 
simulation  results  are  provided  in  Appendix  C. 
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Figure  4.11  -  2D  plume  from  pipeline  break  after  1  year  of  leakage  (VS2DTI). 
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4.2.2    Flare  Pit  Scenario 

The  flare  pit  scenario  was  modelled  for  a  isotropic,  heterogeneous  (coarse-fine-coarse 
layered)  soil  section.  The  flare  pit  was  assumed  to  be  10  m  in  diameter  at  the  base, 
12  m  in  diameter  at  the  surface  and  2  m  deep.  A  constant  head  boundary  was  placed  at 
the  base  of  the  flare  pit  during  periods  when  salt  water  was  assumed  to  be  released  into 
the  pit.  This  scenario  must  therefore  be  modelled  with  an  irregular  (non-rectangular) 
domain  geometry.  The  mesh  generator  provided  with  HYDRUS-2D  only  allows  a 
rectangular  mesh.  An  additional  mesh  generator  is  available  for  purchase,  however,  it 
was  not  purchased  during  this  stage  of  the  study.  For  this  reason,  only  VS2DTI  could  be 
used  to  complete  the  flare  pit  simulation.  Schematics  of  the  scenario  are  shown  in 
Figures  4.11  and  4.12. 


Figure  4.12  -  2D  schematic  of  flare  pit  (red  block  indicates  pit  source). 


Figure  4.13  -  3D  radial  schematic  of  flare  pit  (red  block  indicates  pit  source). 
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The  boundary  conditions  utilized  for  modeling  the  flare  pit  in  an  isotropic,  layered  coarse- 
fine-coarse  soil  section  are  provided  in  Table  4.11.  The  concentration  of  the  water 
discharged  into  the  pit  was  assumed  to  be  equal  to  the  concentration  of  seawater.  The 
salt  water  source  concentration  of  35  kg/m^  was  applied  to  the  base  of  the  pit  for  30  days 
and  then  removed  for  335  days.  This  intermittent  cycle  was  repeated  for  10  years. 


Table  4.11  -  Boundary  conditions  for  flare  pit  scenario. 


Boundary 

Boundary  Type 

Value 

Units 

Surface  Boundary  for  Flow 

Constant  Flux 

5.4  X  10"^ 

m/d 

Radial  Boundary  for  Flow  -  Upper  Aquifer 

Constant  Total  Head 

-3 

m 

Radial  Boundary  for  Flow  -  Lower  Aquifer 

Constant  Total  Head 

-4 

m 

Pit  Boundary  (30  day  discharge  periods) 

Constant  Total  Head 

-1.5 

m 

Pit  Boundary  (335  day  empty  periods) 

Constant  Flux 

5.4x10-^ 

m/d 

Pit  Boundary  for  Solute  (30  day  discharge  periods) 

Solution  Concentration 

35 

kg/m^ 

Pit  Boundary  for  Solute  (335  day  empty  periods) 

Solution  Concentration 

0 

kg/m^ 

The  results  from  the  VS2DTI  model  after  eight  cycles  over  a  period  of  eight  years  are 
provided  in  Figure  4.14.  The  modelling  results  indicates  the  intermittent  recharge  to  the 
pit  has  developed  a  continuous  saturated  zone  in  the  upper  sand  aquifer  and  the  plume 
has  migrated  to  a  radial  distance  of  more  than  50  m.  The  peak  concentration  in  the 
plume  is  close  to  the  applied  source  value  of  35  kg/m^.  The  salt  has  migrated  more  than 
2  m  vertically  into  the  till  beneath  the  surface  aquifer.  Adjacent  to  the  pit,  high  salt 
concentrations  are  predicted  in  the  unsaturated  zone  and  at  the  surface. 
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30  m 


50  m 

Figure  4.14  -  2D  plume  from  a  flare  pit  after  8  years  of  operation  (VS2DTI) 


4.2.3   Salt  Pile  Scenario 


The  salt  pile  was  modelled  for  a  isotropic,  heterogeneous  (coarse-fine-coarse  layered) 
soil  section.  The  salt  pile  was  modelled  as  a  10  m  radius  source.  Salt  was  assumed  to 
enter  the  ground  dissolved  in  surface  infiltration.  Because  of  the  limitations  of  an 
axisymmetric  model,  the  salt  pile  was  limited  to  being  modelled  as  a  circular  surface 
source.  The  source  appears  as  a  disc  on  the  surface  of  the  model,  as  shown  in  Figures 
4.15  and  Figure  4.16. 


Figure  4.15  -  2D  schematic  of  salt  pile  (red  block  indicates  pile). 
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Figure  4.16  -  3D  radial  schematic  of  salt  pile  (red  block  indicates  pile). 

The  boundary  conditions  utilized  for  modelling  the  isotropic,  layered  coarse-fine-coarse 
soil  beneath  the  salt  pile  are  provided  in  Table  4.12.  The  concentration  of  the  infiltrating 
solution  was  assumed  to  be  close  to  seawater  (30  kg/m^).  It  was  assumed  that  leakage 
from  the  salt  pile  continued  at  the  same  location  for  20  years.  After  20  years,  the  source 
was  removed  and  the  scenario  was  simulated  for  an  additional  20  years  allowing  the 
plume  to  dissipate. 


Table  4.12  -  Boundary  conditions  for  salt  pile  scenario. 


Boundary 

Boundary  Type 

Value 

Units 

Surface  Boundary  for  Flow 

Constant  Flux 

5.4  X  10'^ 

m/d 

Radial  Boundary  for  Flow  -  Upper  Aquifer 

Constant  Total  Head 

-3 

m 

Radial  Boundary  for  Flow  -  Lower  Aquifer 

Constant  Total  Head 

-4 

m 

Pile  Boundary  for  Solute  (0  -  20  years) 

Solution  Concentration 

30 

kg/m^ 

Pile  Boundary  for  Solute  (20  -  40  years) 

Solution  Concentration 

0 

kg/m^ 
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The  results  from  the  HYDRUS  and  VS2DTI  models  are  shown  in  Figure  4.17  and 
Figure  4.18.  The  results  indicate  that  the  two  codes  produced  essentially  the  same 
results  for  this  scenario.  Figure  4.17  shows  that  the  salt  plume  has  extended  to  a  depth 
of  approximately  6.0  m,  after  20  years. 
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Figure  4.17  -  Concentration  versus  depth  for  salt  pile  scenario  after  20  years. 


Figure  4.18  (HYDRUS)  and  Figure  4.19  (VS2DTI)  show  that  the  plume  extends  to  a 
radial  distance  of  approximately  5.0  m  away  from  the  site  after  20  years  of  source 
infiltration.  The  contrast  between  the  salt  pile  source  (driven  by  infiltration  alone)  and  the 
"ponded"  flare  pit  source  for  the  same  geology,  shows  that  the  spatial  and  temporal 
distribution  of  the  source  is  critical  information. 
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<  

16.5  m 

Figure  4.19  -  2D  plume  for  salt  pile  after  20  years  of  source  application  (VS2DTI). 
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4.3  Reactive  Transport  Scenarios 

In  many  soils  the  solid  phase  carries  a  net  negative  surface  charge.  For  clay  minerals 
this  charge  is  a  result  of  isomorphous  substitution,  where  structural  cations  of  higher 

valence  are  replaced  by  cations  of  lower  valence  (e.g.  Si4+  is  replaced  by  Al^"*").  As  a 
result,  clay  minerals  have  a  permanent  negative  surface  charge.  Oxides  and  soil 
organic  matter  (SOM)  also  have  charged  surfaces:  With  increasing  pH,  H"*"  is  dissociated 
from  oxide  surfaces  or  from  organic  functional  groups  resulting  in  a  negative  charge. 
Because  dissociation  increases  with  pH  and  ionic  strength,  such  charge  is  variable 
charge  as  opposed  to  the  permanent  charge  of  clays.  At  low  pH  values,  oxides  may 

bind  H"'',  which  results  in  a  positive  surface  charge.  The  pH  values  at  which  positively 
charged  groups  quantitatively  equal  negatively  charged  groups  (i.e.  the  net  surface 
charge  is  zero)  are  called  zero  point  of  charge  (ZPC). 

In  soils,  the  overall  electroneutrality  is  maintained  by  an  excess  of  electrostatically 
attracted  counterions  in  proximity  to  the  charged  surface  (Figure  4.20).  In  the  case  of 
negatively  charged  surfaces,  a  diffuse  double  layer  will  result  where  cations  are  in 
excess  of  anions.  The  excess  ions,  termed  exchangeable  cations,  may  be  exchanged 
with  neutral  salts.  The  quantity  of  exchangeable  cations  (in  meq/kg  dry  soil)  is  defined  as 
the  cation  exchange  capacity  (CEC). 


Distance 


Figure  4.20  -  Distribution  of  ions  near  clay  surface. 
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Exchangeable  cations  are  available  to  plants  (for  example  through  exchange  with  H'' 
liberated  by  the  roots).  Exchange  reactions  are  also  responsible  for  the  retention  of 
introduced  cations  into  the  soil  solution.  In  this  way  the  CEC  gives  the  soil  a  buffering 
capacity,  which  may  slow  down  the  leaching  of  nutrient  cations  and  positively  charged 
pollutants.  CEC  is  usually  dominated  by  the  more  abundant  cations:  Ca,  Mg,  Na,  K,  and 
Al;  thus:  CEC  «  2[Ca]  +  2[Mg]  +  [K]  +  [Na]  +  3[AI]  where  [M]  is  the  meq  per  unit  mass  of 
exchanger  for  exchangeable  ion  M. 

The  most  popular  equations  describing  cation  exchange  are  due  to  Vanselow  (1932), 
Gaines  and  Thomas  (1953),  and  Gapon  (1933).  The  Vanselow  and  Gaines-Thomas 
conventions  differ  only  in  that  Vanselow  is  expressed  in  mole  fractions  and  Gaines- 
Thomas  is  expressed  in  charge  fractions  (equivalents).  The  Gapon  convention,  in 
common  with  the  Gaines-Thomas  convention,  is  also  expressed  in  charge  fractions. 

Because  cation  exchange  is  a  relatively  fast  kinetic  process,  it  can  modify  the  chemical 
composition  of  both  infiltrating  water  and  clays  in  the  solid  phase.  The  exchange 
reaction  of  cations  M  and  N  (with  a  charge  of  m+  and  n+,  respectively),  may  be 
represented  using  the  Gapon  convention  as: 

Adsorbed-Nm  +  nM"^"^  «^  Adsorbed-Mp  +  mN"+ 

Mathematically  this  reaction  has  been  described  by  the  Gapon  equation: 

{ Adsorbed(M)  /  Adsorbed(N) }  =  KnMg  x  { (Mm+)1/m  /  (Nn+)1/n  } 

with  the  left-hand  side  being  the  ratio  of  adsorbed  M  over  N  (both  in  meq  per  mass  unit 
exchanger).  The  right-hand  side  contains  the  reduced  activity  ratio  (where  the  cation 
activities  are  raised  to  a  power  equal  to  the  reciprocal  of  their  valence).  K^m^  in  this 

equation  represents  the  Gapon  selectivity  constant,  which  should  be  constant  over  a 
wide  range  of  conditions. 

The  exchange  reaction  of  cations  M  and  N  (with  a  charge  of  m+  and  n+,  respectively), 
may  also  be  represented  using  the  Gaines-Thomas  convention  as: 


mAdsorbed-N  +  nM"^"*"  o  nAdsorbed-M  +  mN""^ 
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The  Gaines-Thomas  convention  leads  to  the  equation: 

{ Adsorbed(M)n  /  Adsorbed(N)m  }  =  KnMg-t  x  { (M^+)^  I  (Nn+)m  } 

Here  Adsorbed(M)  and  Adsorbed(N)  are  adsorbed  amounts  of  M"^"*"  and  Ni^"*"  ions  in 
meq  per  mass  unit  exchanger.  The  solution  concentrations  on  the  right  hand  side  of  the 
equation  are  the  ion  activities  and  KisiMgyt  '"^  equation  represents  the  Gaines- 
Thomas  selectivity  constant. 

All  three  exchange-equations  (Vanselow,  Gapon,  Gaines-Thomas)  were  developed  for  a 
constant  capacity  (permanent  charge)  exchanger.  This  implies  that  if  the  change  in 
cation  exchange  complex  composition  is  small  compared  to  the  total  store  of 
exchangeable  cations  (which  is  generally  the  case,  particularly  in  the  short  run)  the 
adsorbed  ion  ratio  Adsorbed(N)'"/  Adsorbed(M)"  remains  constant  and  hence  the 
solution  ratio  (N""')"'/(M'"'')"  also  remains  constant. 

The  selectivity  constant  differs  from  unity,  because  small  size  (hydrated)  ions  are 
generally  preferred  over  large  ones  (higher  Knm),  due  to  their  closer  distance  of 
approach  to  the  charged  surface.  The  relative  preferences  (strength  of  retention)  of  the 
most  common  cations,  for  mineral  surfaces,  are  summarized  by  the  Lytropic  Series: 

A|3+  >  Ca2+  >  Mg2+>  K+  >  Na+ 

The  surface  structure  of  the  exchanger  (clay  surface)  may  also  affect  the  selectivity 
constant.  Deviations  from  Lytropic  series  occur  if  specific  chemical  affinity  occurs  for  a 
particular  ion-exchanger  pair. 

UNSATCHEM  uses  the  Gapon  selectivity  constants  for  Ca/Mg,  Ca/Na,  and  Ca/K 
(Simunek  et.  al.  1996).  The  specific  Gapon  constants  are  defined  as  follows: 

KCaMg  =  { Absorbed(Ca)  /  Absorbed(Mg) }  x  {[Mg]"^/  [Ca]"^) 
KCaNa  =  { Absorbed(Ca)  /  Absorbed(Na) }  x  ([Na]  /  [Ca]''') 
KCaK  ={Absorbed(Ca)/Absorbed(K)}  x  ([KJ/ICaf^') 
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where  [M]  denotes  the  equivalent  concentration  in  the  solution  phase  and  Absorbed(M) 
denotes  concentration  in  the  solid  phase. 

VS2DTI  uses  the  Gaines-Thomas  selectivity  coefficients  where  the  constants  for  Ca/Mg, 
Ca/Na,  Ca/K  are  defined  as: 

KCaMg  =  { Absorbed(Ca)2  /  Absorbed(Mg)2  }  x  ([Mg]2  /  [Ca]2) 
KcaNa  =  { Absorbed(Ca)  /  Absorbed(Na)2  }  x  ([Na]2  /  [Ca]) 
KcaK  =  { Absorbed(Ca)  /  Absorbed(K)2  }  x  ([K]2  /  [Ca]) 

Table  4.13  gives  typical  values  for  Gaines-Thonnas  and  Gapon  selectivity  coefficients 
gathered  from  a  brief  review  the  literature. 


Table  4.13  -  Gaines-Thomas  and  Gapon  selectivity  coefficients. 


Selectivity  Coefficient 

Ca-Mg 

Ca-Na 

Ca-K 

Gaines-Thomas 

1.1  -  1.9 

5.0-13 

0.03-0.07 

Gapon 

1.0-1.4 

2.2-3.6 

0.15-0.25 

PHREEQC  and  MINTEQ  Manuals 
Bond  and  Phillips,  1990 
Cernik  etal,  1994 
Schwiech  and  Sardin,  1981 
Valocchi,  1981 
Vulava  etal,  1999 


4.3.1    UNSATCHEM  Cation  Exchange 


Reactive  transport  (cation  exchange)  modelling  was  completed  using  UNSATCHEM  to 
simulate  exchange  reactions  where  solutions  with  high  concentrations  of  Na^  would  be 
infiltrating  into  the  groundwater  system. 

A  ID  column  simulation  was  completed  for  the  coarse-fine-coarse  soil  column  with  an 
initial  groundwater  TDS  composition  of  about  1  kg/m^.  The  infiltrating  solution  had  a 
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concentration  of  30  kg/m^.  Gapon  selectivity  coefficients  used  for  Ca-Mg,  Ca-Na  and 
Ca-K  were  1.2,  2.9  and  0.2,  respectively. 

UNSATCHEM  allows  multicomponent  cation-exchange  models  and  calculates  sodium 
absorption  ratios  (SAR)  defined  as  SAR=  [Na+]  /  { ([Ca2+]  +  [Mg2+])  /  if^.  Figure  4.21 
illustrates  the  sodium  adsorption  ration  (SAR)  profiles  with  time. 


The  initial  SAR  appears  to  be  reasonable  given  the  generic  geology  of  sands  and  tills 
and  the  composition  assumed  for  ambient  groundwater.  When  the  contaminant  solution 
with  a  high  ionic  strength  enters  the  column  with  infiltrating  groundwater,  the  SAR  values 
increase  to  about  500  (consistent  with  the  composition  of  the  salt  water).  UNSATCHEM 
was  designed  to  examine  the  role  of  major  solute  species  in  and  below  the  root  zone  in 
terms  of  irrigation,  fertilization  and  surface  and  groundwater  management.  The  extreme 
salinities  generated  by  salt  wastes  may  be  outside  the  intended  range  of  application  of 
the  model.  Nevertheless,  the  model  does  predict  the  rapid  salinization  of  the  upper 
aquifer  and  can  be  used  to  model  the  response  of  the  system  to  addition  of  Ca-rich 
infiltration  as  a  result  of  surface  application  of  gypsum  as  a  soil  amendment.  Such 
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simulations  may  be  of  greater  value  for  remediation  of  pipeline  breaks  and  flare  pits 
where  Na""  near-surface  concentrations  tend  to  dissipate  over  time. 

4.3.2   VS2DTI  Cation  Exchange 

Reactive  transport  (cation  exchange)  modelling  was  completed  using  VS2DTI  to 
simulate  exchange  reactions  where  solutions  with  high  concentrations  of  Na+  are 
infiltrating  into  the  groundwater  system. 

A  1D  simulation  was  completed  for  the  coarse-fine-coarse  soil  column  with  an  initial 
groundwater  TDS  composition  of  about  1  kg/m^.  The  infiltrating  solution  had  a 
concentration  of  30  kg/m^.  VS2DTI  allows  a  single  ion-pair  cation-exchange  model  to  be 
used  and  saves  only  the  solution  chemistry  as  "total  concentration".  It  was  verified  that 
VS2DTI  was  capable  of  modelling  cation  exchange  scenarios  but  results  were  limited  to 
records  of  bulk  changes  in  mass  balance  for  ions  in  solution  and  on  the  exchanger. 

Any  calculation  of  parameters  such  as  sodium  absorption  ratios,  SAR=  [Na+]  /  { ([Ca2+] 
+  [Mg2+])  /  2}^'^,  would  require  source  code  modifications  to  VS2DTI  to  write  additional 
information  to  output  files.  Some  kind  of  script  could  then  be  written  to  post  process  the 
additional  VS2DTI  output  files.  In  this  regard,  UNSATCHEM  is  significantly  more  useful 
for  cation-exchange  modelling  in  the  present  salt-contamination  context. 

Because  VS2DTI  has  2D  Tier  3  application  potential,  the  cation  exchange  capabilities 
may  eventually  prove  to  be  more  valuable  than  the  simple  1-D  column  model  of 
UNSATCHEM.  At  the  present  stage  UNSATCHEM  is  recommended. 
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5.0  RECOMMENDATIONS 

5.1  Responsible  Numerical  Modelling 

Numerical  models  can  be  useful  tools  for  understanding  contaminant  transport  at  a  site. 
However,  the  results  from  the  numerical  model  are  dependent  upon  factors  such  as 
input  data,  mesh  design,  convergence  criteria,  and  mass  balance. 

5.1.1  Input  Data 

Input  data  should  be  based  primarily  on  site-specific  data.  Where  site-specific  data 
cannot  be  obtained,  it  may  be  adequate  to  input  data  based  on  parameter  estimates 
from  the  literature.  The  adequacy  of  the  geologic  model,  material  properties, 
contaminant  transport  properties,  and  groundwater  flow  system  can  be  evaluated  by 
viewing  the  model  calibration. 

Calibration  targets  can  be  pre-determined  based  on  generally  accepted  criteria.  Both 
qualitative  and  quantitative  calibration  criteria  can  be  utilized  in  the  calibration  process. 

The  following  list  constitutes  a  set  of  qualitative  targets  that  an  "acceptable"  model 
should  satisfy: 

•  Plume  shape  and  extent  estimates  based  on  concentration  contours; 

•  Observed  downward  movement  of  brine  front  using  concentration  contours; 

•  Observed  hydraulic  gradients  and  flow  directions; 

•  Spatial  distribution  of  groundwater  highs  and  lows;  and, 

•  Estimated  fluid  fluxes  in  all  aquifers. 

5.1.2  Mesh  Design 

A  poorly  designed  mesh  can  result  in  both  numerical  dispersion  and  numerical 
oscillation.  One  indication  of  numerical  problems  is  the  development  of  unreasonable 
concentrations  (such  as  concentrations  exceeding  source  concentration).  Numerical 
oscillation  and  numerical  dispersion  can  be  controlled  by  satisfying  mesh  size  and  time 
step  constraints  imposed  by  the  Peclet  and  Courant  number  criteria. 
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Satisfying  the  Peclet  number  criterion  ensures  that  the  elements  are  small  enough  to 
avoid  dispersion  that  occurs  because  the  movement  of  contaminant  from  cell-to-cell  in 
the  mesh  as  a  result  of  successive  steps  in  time  occurs  faster  than  the  physical  mixing 
process.  Numerical  dispersion  is  an  artefact  of  the  mesh.  The  smallest  distance  that 
the  contaminant  can  travel  in  one  time-step  is  the  distance  between  adjacent  mesh 
nodes;  if  the  time-step  or  mesh  spacing  is  too  large,  numerical  dispersion  can  become 
the  dominant  mixing  process.  A  high  degree  of  spatial  resolution  is  required  to  avoid 
excessive  "artificial"  mixing.  The  Peclet  constraint  for  the  grid  size  requires  that: 


where  Dl  and  ai  are  respectively  the  longitudinal  dispersion  coefficient  and  longitudinal 
dispersivity.  This  requires  that  the  mesh  size  in  the  direction  of  flow  is  less  than  2  x  the 
dispersivity. 

Satisfying  the  Courant  number  criterion  ensures  that  the  distance  of  advective 
contaminant  transport  in  a  single  time  step  does  not  exceed  one  grid  block  or  element 
length.  The  Courant  constraint  for  the  time-step  size  requires  that: 


where  v  is  the  advective  velocity.  Remembering  v  =  -  K(Ah/Ax)/n  where  K  is  hydraulic 
conductivity,  n  is  porosity  and  Ah/Ax  is  the  hydraulic  gradient.  It  is  easy  to  see  that  At 
actually  has  a  quadratic  rather  than  linear  dependence  on  Ax.  Unfortunately,  when  a 
fine  grid  is  used  to  reduce  numerical  dispersion,  the  time  step  for  dispersive  transport 
calculations  may  become  smaller  than  the  ideal  time  step  for  advective  calculations, 
because  the  dispersive  step  to  satisfy  the  Courant  constraint  has  quadratic  dependence 
on  grid  size. 

Numerical  instabilities  (oscillations)  in  the  calculation  of  diffusion/dispersion  are 
eliminated  with  the  Von  Neumann  condition  constraint  (really  a  combination  of  the  Peclet 
and  Courant  constraints): 


Ax  <  2Dl  /  V  «  2aL 


At<Ax/  V 


At  <  Ax^  /  3Dl  =  Ax^  /  SqlV 
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It  is  easy  to  see  that  at  low  advective  velocities  the  constraint  is  readily  satisfied  but  that 
progressively  smaller  At  is  required  as  Ax  decreases. 

Spatial  and  temporal  discretization  of  the  model  domain  must  be  designed  carefully  to 
represent  a  natural  system  and  avoid  artifacts  created  by  the  numerical  dispersion  and 
oscillations  in  the  calculations.  Artificial  "mixing"  or  numerical  dispersion  is  the  most 
serious  problem  encountered  in  dispersive/diffusive  transport  modelling  of  dense  brine 
transport  problems  (Uwiera,  1998). 

Appendix  D  illustrates  the  result  of  a  VS2DTI  solution  where  numerical  oscillation  has 
occurred. 

5.1.3    Convergence  Criteria  and  IVIass  Balance 

Problems  with  the  specification  of  convergence  criteria  can  often  be  diagnosed  if  the 
mass  balance  of  the  model  exceeds  an  acceptable  tolerance.  Reducing  the 
convergence  criteria  for  head  and/or  concentration  can  usually  decrease  the  mass 
balance  error  in  a  particular  simulation.  Other  factors  that  may  affect  the  mass  balance 
are  mesh  spacing  and  time  step  size.  Reducing  these  values  should  also  lead  to  an 
improved  mass  balance,  provided  the  Peclet,  Courant  and  Von  Neumann  numerical 
constraints  continue  to  be  honoured. 

There  are  situations  when  a  relatively  large  mass  balance  error  may  not  be  indicative  of 
solution  error.  Changing  constant  head  boundary  or  initial  conditions  to  simulate 
changes  in  site  configuration  may  produce  large  mass  balance  errors.  The  model  sees 
this  as  an  instantaneous  change  in  head  and  accompanying  this  head  change  will  be  an 
instantaneous  change  in  fluid  stored  within  an  element.  Since  this  increase  (or 
decrease)  in  storage  was  not  accounted  for  in  the  previous  time  step  the  model  treats  it 
as  a  discrepancy  in  the  total  mass  balance. 

A  common  occurrence  of  this  apparent  discrepancy  is  on  the  very  first  time  step  of  a 
simulation  when  there  is  an  inconsistency  between  change  in  initial  and/or  boundary 
conditions.  Such  an  occurrence  is  manifested  by  a  large  mass  balance  error  for  the  first 
time  step. 
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For  this  reason  it  is  difficult  to  recommend  a  mass  balance  criterion  that  all  simulations 
must  satisfy.  In  the  absence  inconsistencies  in  initial  conditions  and  boundary  conditions 
or  abrupt  changes  in  boundary  conditions,  a  mass  balance  error  of  less  than  5%  is  a 
reasonable  target. 

5.2  Discussion  of  Recommended  Codes 

Four  of  the  codes  (CHEMFLO,  HYDRUS,  UNSATCHEM,  and  VS2DTI)  can  be  readily 
applied  to  1D  Tier  2  analysis  for  the  fate/transport  of  salt.  All  four  codes  produced 
comparable  consistent  results  for  the  Tier  2  generic  scenarios. 

SEVIEW  (AT  123D,  BIOSCREEN  and  SESOIL)  was  the  only  code  tested  that  was 
designed  to  predict  groundwater  recharge  from  seasonal  climatic  data.  SESOIL  is  not 
recommended  because  of  the  difficulties  encountered  in  obtaining  reasonable  estimates 
of  groundwater  recharge  through  tills  for  the  arid  western  prairie  environment. 
Nevertheless,  AT123D  can  be  used  to  carry  out  1D  Tier  2  analysis  but  without  SESOIL 
has  no  advantages  over  the  other  codes  and  is  limited  to  homogeneous  (or  equivalent 
homogeneous)  materials. 

For  the  more  complex  Tier  3  analysis  only  VS2DTI  and  HYDRUS-2D  were  capable  of 
handling  axisymmetric  problems.  Both  codes  performed  consistently  on  the  2D  generic 
scenarios  but  VS2DTI  proved  more  flexible  in  use. 

For  cation  exchange  problems  only  UNSATCHEM  and  VS2DTI  provided  a  means  of 
calculation.  The  two  codes  use  different  definitions  for  selectivity  coefficients  and  are 
difficult  to  compare  head-to-head.  UNSATCHEM  is  the  most  sophisticated  reactive 
transport  model  and  readily  predicts  soil  quality  parameters  such  as  sodium  absorption 
ratio  (SAR)  and  electrical  conductivity  (EC). 

5.2  CHEMFLO 

CHEMFLO  is  likely  the  simplest  of  the  codes  tested.  CHEMFLO  is  a  public  domain  code 
that  is  supported  by  Oklahoma  State  University  (OSU)  Oklahoma  Agricultural 
Experiment  Station.  The  version  of  CHEMFLO  used  for  the  tests  was  CHEMFLO-2000 
beta  version  2003.02.26  supplied  by  the  author  with  added  functionality  prior  to  public 
release. 
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CHEMFLO  is  an  excellent  and  effective  tool  for  Tier  2  analysis.  General  knowledge  of 
the  stratigraphic  section,  elevations  of  potentiometric  surfaces,  groundwater  recharge 
rate,  basic  material  properties,  and  basic  chemical  transport  properties  are  the  minimum 
required  parameters  to  run  CHEMFLO. 

CHEMFLO's  main  strength  is  its  ability  to  provide  results  for  a  simple  Tier  2  analysis 
relatively  quickly  with  simple  data  input  requirements.  The  input  parameters  are 
relatively  straightfonA/ard  and  well  documented  and  the  user  interface  is  intuitive  and 
easy  to  use.  CHEMFLO  also  has  a  java-based  interface  that  gives  it  added  flexibility  in 
terms  of  platform  independence. 

CHEMFLO's  weakness  is  that  it  is  a  1D  column  model.  The  current  public  release 
(2002.08.08)  of  CHEMFLO  can  only  model  the  saturated  zone  to  1  m  below  the  water 
table.  A  second  problem  with  the  current  release  of  CHEMFLO  (2002.08.08)  is  the 
restrictions  it  places  on  the  Van  Genuchten  beta  parameter.  Both  these  restrictions  are 
removed  in  the  version  tested  in  this  study  and  the  restrictions  will  not  be  present  in  the 
next  public  release  of  the  code.  CHEMFLO  is  also  limited  to  simulating  a  maximum  of 
four  materials,  however,  this  should  be  adequate  enough  for  most  simple  Tier  2  analysis 
problems. 

An  annoyance  associated  with  CHEMFLO  is  that  it  does  not  save  the  input  data 
between  runs.  Printing  the  screen  to  the  clipboard  and  saving  it  in  a  word  processing 
program  such  as  MS-Word  can  solve  this  problem  relatively  easily,  but  the  data  must  be 
re-entered  every  time  the  program  is  executed  rather  than  read  from  a  file. 

Although  the  current  release  of  CHEMFLO  (2002.08.08)  has  limitations,  the  author 
provided  a  beta  version  of  CHEMFLO  (2003.02.26)  to  MDH  that  corrected  the  limitation 
of  the  Van  Genuchten  beta  parameter  and  allowed  for  a  larger  saturated  zone  to  be 
modelled.  The  author  expected  that  the  new  public  release  (complete  with 
documentation)  will  be  available  in  March,  2003,  however,  no  definite  time  was  given  for 
the  release  date. 

CHEMFLO  was  verified  against  analytical  solutions  and  other  numerical  codes  as  part  of 
the  evaluation  process  and  appears  to  produce  accurate  consistent  results.  The  code  is 
reliable,  well  documented,  easy  to  use,  in  the  public  domain  and  supported  by  a 
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reputable  technical  institution.  With  a  java-based  interface,  it  is  also  has  the  added 
advantage  of  platform-independence. 

The  revised  version  of  CHEI\/IFLO  to  be  released  after  l\/larch  2003  is  recommended 
as  first  choice  for  Tier  2  analysis. 

5.3  HYDRUS 

HYDRUS-2D  was  one  of  two  2D  codes  tested  in  the  group  of  five  codes  selected  for 
detailed  review  in  Phase  II.  HYDRUS  is  a  proprietary  code  that  is  supported  by  the 
United  States  Salinity  Laboratory  (USSL)  and  distributed  by  the  International 
Groundwater  Modelling  Centre  (IGWMC)  at  the  Colorado  School  of  Mines.  The  version 
tested  was  HYDRUS-2D  version  2.007. 

HYDRUS  is  an  effective  tool  for  Tier  2  analysis  and  can  be  applied  to  most  Tier  3 
problems.  General  knowledge  of  the  stratigraphic  section,  elevations  of  potentiometric 
surfaces,  groundwater  recharge  rate,  basic  material  properties,  and  basic  chemical 
transport  properties  are  the  minimum  required  parameters  to  run  HYDRUS  for  1D 
column  analysis  and  for  2D  axisymmetric  problems. 

The  main  strength  of  HYDRUS  is  its  ability  to  provide  results  for  both  simple  Tier  2 
analysis  and  more  complex  2D  axisymmetric  analysis  for  Tier  3  assessment.  HYDRUS 
automatically  adjusts  the  time  step  to  ensure  that  the  Courant  criterion  is  satisfied  and 
reports  the  critical  grid  Peclet  number.  This  is  a  valuable  feature  for  ensuring  the  quality 
of  results. 

The  main  disadvantages  of  HYDRUS  are  its  inability  to  simulate  an  irregular  (non- 
rectangular)  geometry  without  purchasing  an  additional  proprietary  mesh  generator  and 
a  large  number  of  "built-in"  undocumented  restrictions. 

1.  With  the  "basic"  mesh  generator  supplied  with  HYDRUS,  specifying  a  mesh  with 
variable  spacing  is  tedious.  When  the  mesh  is  changed  material  property 
boundaries  and  boundary  conditions  must  be  re-specified.  The  version  of 
HYDRUS  tested  (2.007)  crashes  if  the  mesh  spacing  is  expanded  using  an 
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expansion  factor  above  a  critical  value  (but  does  not  warn  the  user  of  such  a 
restriction). 

2.  The  HYDRUS  solver  stability  also  appears  to  be  somewhat  sensitive  to  mesh 
spacing.  Again  the  program  crashes  without  prior  warning  of  a  restriction  and  the 
user  may  need  to  increase  or  decrease  the  mesh  size  (and  must  consequently  re- 
specify  all  material  boundaries  and  flow  and  transport  boundary  conditions). 

3.  For  the  generic  axisymmetric  scenario  with  a  pipe  break,  it  was  necessary  to  use 
an  "internal"  nodal  flux  boundary  condition  with  a  specified  source  concentration. 
Such  boundary  conditions  are  poorty  documented  in  HYDRUS  and  there  is  an 
unspecified  maximum  flux  rate  that  can  be  applied.  Again  the  program  crashes 
without  warning  if  this  threshold  value  is  exceeded. 

4.  HYDRUS  cannot  simulate  multiple  recharge  periods  with  different  boundary 
conditions  within  the  same  simulation.  However,  the  package  does  allow  the  user 
to  import  concentrations  and  pressure  heads  from  previous  simulations. 
Unfortunately,  there  is  an  unspecified  restriction  on  the  number  of  nodal 
concentration  values  that  can  be  imported  in  an  ASCII  file  (although  no  such 
restriction  applies  to  heads  or  to  concentrations  if  the  file  is  saved  in  binary  format). 

HYDRUS-2D  is  certainly  adequate  to  use  for  Tier  2  analysis  but  has  little  to  offer  over 
CHEMFLO,  UNSATCHEM  and  VS2DTI  in  this  role.  HYDRUS-1D  can  also  perform  the 
same  function  (but  was  not  tested  in  this  phase  of  the  study).  Over  all,  HYDRUS  is  a 
well-written  and  reliable  numerical  code  that  produced  accurate  results  on  all  verification 
problems.  It's  major  fault  is  a  poor  user  interface  that  fails  to  protect  the  user  from 
predictable  parameter  limitations  and  input-output  restrictions. 

Although  the  HYDRUS-2D  code  is  fundamentally  sound  and  has  many  excellent 
features,  the  limitations  of  the  program  are  significant  enough  to  not 
recommended  its  use  for  Tier  3  analysis. 

5.4  SEVIEW 

SEVIEW  is  a  proprietary  "GUI  wrapper"  around  the  public-domain  BIOSCREEN, 
SESOIL  and  AT123D  codes.   The  GUI  generates  input  for  the  AT123D  contaminant 
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transport  code  using  SESOIL.  The  version  of  the  package  tested  was  SEVIEW 
version  6.0. 

Detailed  knowledge  of  the  seasonal  climatic  data  and  chemical  properties  are  required 
to  run  SESOIL.  General  knowledge  of  the  stratigraphic  section,  groundwater  flow 
conditions,  groundwater  recharge  rate  and  basic  soil  properties  are  the  minimum 
required  parameters  to  run  Tier  2  1D  column  problems  in  SEVIEW. 

SEVIEW's  main  strength  is  the  ability  of  the  SESOIL  model  to  predict  groundwater 
recharge  from  measured  monthly  or  daily  climatic  data.  The  performance  of  SESOIL  as 
a  tool  to  generate  groundwater  recharge  from  a  water  balance  was  disappointing. 

SEVIEW's  main  weakness  is  SESOIL's  requirement  of  extensive  climatic  and  chemical 
input  data.  Using  "raw"  monthly  climate  data  for  sites  in  the  western  Canadian  prairies 
SESOIL  was  not  able  to  generate  credible  groundwater  recharge  rates  for  fine-grained, 
low-permeability  soils  (despite  numerous  attempts).  SEVIEW  allows  up  to  four 
heterogeneous  soil-layers  but  generates  a  single  "equivalent"  homogeneous  material  for 
the  SESOIL  analysis.  Only  a  homogeneous  material  can  be  input  into  AT123D.  The 
inability  of  SEVIEW  to  simulate  truly  heterogeneous  soil  systems  limits  the  number  of 
sites  in  Alberta  where  the  code  could  be  applied  (even  for  Tier  2  analysis).  SEVIEW 
appears  to  be  designed  for  organic  contaminants  at  low  concentrations  in  moderate  to 
high  permeability  agricultural  soils. 

Without  SESOIL,  AT123D  has  many  limitations  and  no  advantages  over  CHEMFLO, 
HYDRUS,  UNSATCHEM  and  VS2DTI.  For  this  reason,  systematic  further  testing  of 
AT123D  was  discontinued. 

SEVIEW  is  not  recommended  as  a  suitable  code  for  screening  of  salt 
contaminated  sites  in  Alberta. 

5.5  UNSATCHEM 

UNSATCHEM  is  a  public  domain  code  that  was  written  by  the  United  States  Salinity 
Laboratory  and  is  receiving  limited  support.  It  has  been  superseded  by  HYDRUS-1D 
(which  was  not  tested).  The  version  of  UNSATCHEM  used  for  the  tests  was 
UNSATCHEM  version  2.0  with  a  compilation  date  in  1996. 
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General  knowledge  of  the  stratigraphic  section,  elevations  of  potentiometric  surfaces, 
groundwater  recharge  rate,  basic  material  properties,  and  basic  chemical  transport 
properties  are  the  minimum  required  parameters  to  run  UNSATCHEM. 

UNSATCHEM  is  a  reliable  and  relatively  easy-to-use  1D  code  for  Tier  2  analysis  but  is 
not  the  first-choice  for  such  applications. 

UNSATCHEM's  main  advantage  lies  in  it  ability  to  simulate  cation  ion-exchange 
reactions  as  a  component  of  simple  Tier  2  1D  column  problems.  The  input  parameters 
are  relatively  straightforward  and  reasonably  well  documented  (although  some  sign 
conventions  are  poorly  explained).  The  user  interface  is  less  than  intuitive  but  relatively 
easy-to-use  once  familiarity  has  been  gained. 

UNSATCHEM's  main  disadvantages  are  associated  with  the  specification  of  boundary 
conditions.  A  constant  flux  boundary  must  be  applied  as  negative  for  movement  of 
water  into  the  column  and  positive  for  movement  of  water  out  of  the  column.  This  flux 
boundary  condition  is  recommended  for  use  over  the  atmospheric  flux  boundary 
condition  because  UNSATCHEM  appears  to  treat  the  concentration  flux  boundary  as  a 
constant  concentration  boundary  when  an  atmospheric  flux  boundary  condition  is 
applied  at  the  surface. 

Another  disadvantage  of  UNSATCHEM  is  its  inability  to  import  pressure  head  or 
concentration  data  from  previous  runs  or  to  use  multiple  recharge  periods  with  a 
constant  flux  boundary.  In  order  to  complete  the  1D  simulations,  MDH  had  to  manually 
input  the  output  data  from  one  run  as  the  initial  condition  for  a  second  run  to  simulate  the 
removal  of  a  source.  This  inflexibility  limits  the  usefulness  of  UNSATCHEM. 

A  major  concern  with  UNSATCHEM  is  the  lack  of  further  development  and  the  limited 
nature  of  the  support  for  the  code.  HYDRUS-1D  is  the  proprietary  replacement  for  the 
public  domain  UNSATCHEM  code  but  lacks  the  sophistication  of  UNSATCHEM  for 
reactive  transport  modelling. 

In  order  to  use  the  advanced  reactive  transport  capabilities  of  UNSATCHEM,  addition 
information  on  the  chemical  parameters  (particularly  soil  chemistry  and  mineralogy, 
cation  exchange  capacity,  and  ambient  groundwater  and  pore-water  chemistry)  is 
required.  Such  data  is  generally  collected  as  part  of  Tier  3  assessments.  The  reactive 
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transport  capabilities  of  UNSATCHEM  can  provide  a  means  of  quantifying  the  effects  of 
applying  gypsum  amendment  to  the  surface  soils  in  the  vicinity  of  remediated  pipe 
breaks  and  flare  pit  sites. 

UNSATCHEM  was  verified  against  analytical  solutions  and  other  numerical  codes  as 
part  of  the  evaluation  process  and  appears  to  produce  accurate  consistent  results.  The 
code  is  reliable,  well  documented,  easy-to-use,  and  in  the  public  domain. 

UNSATCHEM  is  recommended  for  use  in  Tier  3  assessment  wiiere  simulation  of 
cation  exchange  processes  is  a  requirement 

5.6  VS2DTI 

VS2DTI  is  the  second  of  two  2D  codes  tested  in  Phase  2  of  the  project.  VS2DTI  is  a 
public  domain  code  supported  by  the  United  States  Geological  Survey  (USGS).  The 
version  tested  was  VS2DTI  version  1.1. 

VS2DTI  is  an  excellent  code  for  both  Tier  2  (ID  columns)  and  Tier  3  (2D  axisymmetric) 
analysis  of  salt  transport  (at  low  concentration).  General  knowledge  of  the  stratigraphic 
section,  elevations  of  potentiometric  surfaces,  groundwater  recharge  rate,  basic  material 
properties,  and  basic  chemical  transport  properties  are  the  minimum  required 
parameters  to  run  VS2DTI  for  1 D  column  analysis  and  for  2D  axisymmetric  problems. 

VS2DTrs  main  advantage  lies  in  its  extreme  flexibility  (it  was  the  only  code  to  simulate 
every  one  of  the  generic  scenarios  without  significant  difficulty).  VS2DTI  has  a  good 
user-friendly  interface  although  some  of  the  important  time  saving  features  are  not 
obvious  to  the  first-time  user.  VS2DTI  separates  the  definition  of  the  model  domain, 
material  properties  and  boundary  conditions  from  the  generation  of  a  mesh.  This  means 
that  a  problem  can  be  re-meshed  in  seconds  and  identical  problems  are  readily  run  on 
multiple  meshes.  This  is  a  major  advantage  when  checking  the  validity  of  results. 
VS2DTI  also  allows  multiple  recharge  periods  to  be  specified  with  different  boundary 
conditions  and  numerical  time-stepping  parameters.  Again,  this  provides  for  easy 
checking  and  adds  to  flexibility.  VS2DTI  has  the  ability  to  model  cation  exchange 
reactions  although  the  post-processing  options  in  respect  of  this  feature  are  minimal. 
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VS2DTrs  main  disadvantage  is  that  solutions  can  be  extremely  sensitive  to  mesh 
discretization  and  sometimes,  initial  conditions.  Unlike  HYDRUS,  VS2DTI  does  not 
automatically  adjust  the  time  step  to  satisfy  the  Courant  condition  and  does  not 
automatically  output  the  grid  Peclet  number.  Extreme  care  is  needed  with  mesh  design 
and  the  specification  of  initial  conditions  to  ensure  that  "valid"  results  are  obtained. 
VS2DTI  does  report  mass  balance  errors  to  the  screen  during  simulation  runs  but  their 
interpretation  requires  technical  knowledge  and  experience. 

Although  VS2DTI  allows  the  user  to  input  the  number  of  iterations  for  convergence  of 
transport  steps,  the  program  ignores  the  parameters  input  by  the  user  if  a  threshold 
value  of  200  is  exceeded.  Editing  and  re-compiling  the  source  code  could  easily  fix  this, 
but  the  interface  should  warn  the  user  of  the  limitation.  Another  disadvantage  of  VS2DTI 
is  the  inability  to  cut  cross-sections  of  data  and  to  view  graphs  of  concentration  or  head 
versus  time  or  distance  in  the  postprocessor.  The  CHEMFLO,  HYDRUS  and 
UNSATCHEM  post-processing  capabilities  are  better  in  this  area. 

VS2DTI  models  cation  exchange  using  the  Gaines-Thomas  convention  for  selectivity 
coefficients.  Unfortunately,  the  program  appears  to  save  only  rudimentary  information 
about  reactive  chemical  processes.  It  is  not  possible  to  plot  cation  ratios  either  against 
time  or  spatially  for  either  the  solution  or  the  exchanger.  Source  code  modifications 
would  be  necessary  to  implement  such  a  feature. 

VS2DTI  was  verified  against  analytical  solutions  and  other  numerical  codes  as  part  of 
the  evaluation  process  and  appears  to  produce  accurate  consistent  results.  The  code  is 
reliable,  well  documented,  easy  to  use,  in  the  public  domain  and  supported  by  a 
reputable  technical  institution.  With  a  java-based  interface,  it  is  also  has  the  added 
advantage  of  platform-independence 

VS2DTI  is  recommended  for  application  to  boti)  1D  column  (Tier  2)  and  2D 
axisymmetric  (Tier  3)  scenarios. 
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6.0  DISCLAIMER 

MDH  Engineered  Solutions  Corp.,  hereinafter  collectively  referred  to  as  "MDH",  has 
exercised  reasonable  skill,  care  and  diligence  in  assessing  the  modelling  codes  within 
this  report,  but  makes  no  guarantees  or  warranties  as  to  the  accuracy  or  completeness 
of  this  assessment.  While  all  reasonable  efforts  have  been  made  to  verify  the  accuracy 
of  the  information  upon  which  this  evaluation  is  based,  MDH  Engineered  Solutions  Corp. 
will  not  be  liable  under  any  circumstances  for  the  direct  or  indirect  damages  incurred  by 
any  individual  or  entity  due  to  the  contents  of  this  report,  omissions  and/or  errors  within, 
or  use  thereof,  including  damages  resulting  from  loss  of  data,  loss  of  profits,  loss  of  use, 
interruption  of  business,  indirect,  special,  incidental  or  consequential  damages,  even  if 
advised  of  the  possibility  of  such  damage.  This  limitation  of  liability  will  apply  regardless 
of  the  form  of  action,  whether  in  contract  or  tort,  including  negligence. 

MDH  Engineered  Solutions  Corp.  has  prepared  this  report  for  the  exclusive  use  of 
Alberta  Environment  and  does  not  accept  any  responsibility  for  the  use  of  this  report  for 
any  purpose  other  than  that  intended.  Any  alternative  use,  reliance  on,  or  decisions 
made  based  on  this  document  are  the  responsibility  of  the  alternative  user  or  third  party. 
MDH  Engineered  Solutions  Corp.  accepts  no  responsibility  to  any  third  party  for  the 
whole  or  part  of  the  contents  and  exercises  no  duty  of  care  in  relation  to  this  report. 
MDH  has  assessed  each  modelling  code  for  the  specific  purpose  required  by  Alberta 
Environment  and  accepts  no  responsibility  for  damages  suffered  by  any  third  party  as  a 
result  of  decisions  made  or  actions  based  on  this  report. 

The  mention  of  a  tradename  is  solely  for  illustrative  purposes.  MDH  does  not  hereby 
endorse  any  tradename,  warrant  that  a  tradename  is  registered,  or  approve  a  tradename 
to  the  exclusion  of  other  tradenames.  MDH  does  not  give,  nor  does  it  imply,  permission 
or  license  for  the  use  of  any  tradename. 
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7.0  CLOSURE 

We  trust  that  this  report  meets  all  your  present  requirements.  Should  you  have  any 
questions  or  comments  please  contact  us.  We  look  forward  to  discussing  this  report 
further  with  you  in  the  near  future. 

Respectfully  submitted, 
MDH  Engineered  Solutions 


Association  of  Professional  Engineers, 
Geologists,  and  Geophysists  of  Alberta 
Permit  to  Practice  P7607 


Roxanne  Pauls,  M.Sc. 


Andrew  Karvonen,  M.Sc,  P. Eng.,  P.Geo.  (Sask) 


Dr.  Malcolm  Reeves,  P. Eng.,  P.Geo.  (Sask) 


Dr.  Moir.D.  Haug,  P. Eng. 


^MDH 
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Table  A2  Complete  input  data  requirements  for  CHEMFLO,  HYDRUS,  UNSATCHEM,  VS2DTI 


This  table  identifies  what  parameters  need  to  be  input  with  each  model.  If  the  cell  says  Yes"  then  the  variable  must  be  entered  into  the 
model.  If  the  cell  says  "No"  then  there  is  no  option  to  enter  the  variable  into  the  model.  If  the  cells  sayOpf  then  the  variable  can  be  input 
into  the  model,  but  nothing  needs  to  be  input  into  the  model  if  it  is  not  clicked  on  in  the  model. 


I Must  be  measured/providec 
Must  be  estimated 
Not  needed/Defaultee 


# 

Variable 

Unit 

Program 

CHEMFLO 

SEVIEW 

UNSATCHEM 

HYDRUS2D 

VS2DTI 

GEOLOGY: 

"T" 

^^HKabove  bedrock  surface 

Yes 

Yes 

Yes 

Yes 

Yes 

2 

^^^B  below  bedrock  surface 

Yes 

Yes 

Yes 

Yes 

Yes 

3 

^^Hibelow  zone  of  contamlnatlor 

Yes 

Yes 

Yes 

Yes 

Yes 

SOIL  PROPERTIES: 

^H^^^HHte  content 

- 

Yes 

No 

Yes 

Yes 

Yes 

5 

HPHPIPIIPto' Content 

- 

Yes 

No 

Yes 

Yes 

Yes 

6 

- 

No 

Yes 

No 

No 

No 

7 

Permeability 

Saturated  hydraulic  conductivity 

L/T 

Yes 

Yes 

Yes 

Yes 

Yes 

 .i....                            Intrinsic  permeability 

L2 

No 

Yes 

No 

No 

No 

~8~ 

M/L3 

Yes 

Yes 

Yes 

Yes 

Yes 

SOLUTE  PROPERTIES: 

T 

10 

I^^R^l^K  capacity 

mol/M 

No 
No 

Yes 
Yes 

No 
Yes 

No 
No 

No 

FIELD  CONDITIONS 

Opt 

"T7 

HHbiessure  head  at  column  surface 

L 

Yes 

No 

Yes 

Yes 

Yes 

12 

nBl^ressure  head  data  above  water  table 

L 

Yes 

No 

Yes 

Yes 

Yes 

13 

EiSDcitlon  of  water  table 

L 

Yes 

No 

Yes 

Yes 

Yes 

Initial  pressure  head  data  below  water  table 

L 

No 

No 

Yes 

Yes 

Yes 

15 

Pressure  head  in  all  aquifers 

L 

Yes 

No 

Yes 

Yes 

Yes 

16 

Location  of  contamination  source: 

Yes 

Yes 

Yes 

Yes 

Yes 

17 

Water  analysis 

Yes 

Yes 

Yes 

Yes 

Yes 

18 

Soli  salinity  anajysis  with  depth 

Yes 

Yes 

Yes 

Yes 

Yes 

19 

initial  concentration  at  column  surface 

M/L3 

Yes 

Yes 

Yes 

Yes 

Yes 

# 

Unit 

Program 

Additional  2-D  parametei 

CHEMFLO 

SEVIEW 

UNSATCHEJVI 

HYDRUS2D 

VS2DTI 

20 

Detailed  site  plar  ^ 

No 

Yes 

No 

Yes 

Yes 

21 

Topographic  map  of  are« 

No 

Yes 

No 

Yes 

Yes 

22 

Indication  of  receptors  away  from  site  boundarie 

No 

Yes 

No 

Yes 

Yes 

23 

Conceptual  Idea  of  lateral  variation  in  geolog 

No 

Yes 

No 

Yes 

Yes 

24 

Conceptual  idea  of  lateral  variation  of  water  tabi 

No 

Yes 

No 

Yes 

Yes 

25 

Conceptual  Idea  of  lateral  variation  in  pressure  hea 

No 

Yes 

No 

Yes 

Yes 

# 

Variable 

Unit 

Program 

CHEMFLO 

SEVIEW 

UNSATCHEM 

HYDRUS2D 

VS2DTI 

SOIL  PROPERTIES: 

26 

SWCC 

- 

Qm  (soil  water  retention  function 

- 

No 

No 

Yes 

Opt 

No 

Qa  {soil  water  retention  function 

No 

No 

Yes 

Opt 

No 

Alpha  parameter 

1/L 

Yes 

No 

Yes 

Yes 

Yes 

Beta  or  n  parametei 

Yes 

No 

Yes 

Yes 

Yes 

27 

Soil  pore  disconnectness 

No 

Yes 

No 

No 

No 

28 

Specific  storage 

No 

No 

No 

No 

Yes 

SOLLUTE  PROPERTIES: 

29 

Dispersion  coefflcleni 

L2n- 

30 

longitudinal  dispersivity 

L 

Yes 

Yes 

Yes 

Yes 

Yes 

31 

Transverse  dispersivity 

L 

No 

Yes 

No 

Yes 

Yes 

32 

Diffusion  coefficient  (ions  in  porous  medum 

L2/T 

Yes 

Yes 

Yes 

Yes 

Yes 

33 

Ko 

Uniform  soil  water  partition  coefficien 

L3/M 

Yes 

Yes 

No 

No 

kd 

Freundlich  adsorption  isothem 

No 

Yes 

No 

No 

off 

Adsorption  isotherm  coefficient  k; 

No 

No 

No 

Yes 

No 

34 

Ion  Exchange 

Km:  Gaines-Thomas  ion-exchange  selectivity  coefficient 

No 

No 

No 

No 

Opt 

Kg:  Gapon  ion-exchange  selectivity  coefficient 

No 

No 

Yes 

No 

No 

35 

Solubility  In  water 

M/L3 

No 

Yes 

No 

No 

No 

36 

Molecular  weight 

M/mol 

No 

Yes 

No 

No 

No 

37 

IHenry's  constant 

No 

Yes 

No 

Yes 

Opt 

38 

Chemical  valence 

No 

Yes 

No 

No 

No 

BOUNDARY  CONDITIONS 

39 

Hydraulic  Surface  Boundary 

Constant  flux 

UT 

Opt 

No 

Opt 

Opt 

Opt 

Variable  flux 

LfT 

No 

Opt 

Opt 

Opt 

Opt 

Constant  heac 

L 

Opt 

No 

Opt 

Opt 

Opt 

Atmospheric  BC  with  surface  layei 

LTT 

No 

No 

Opt 

Opt 

No 

Atmospheric  BC  with  surface  runof 

LfT 

No 

No 

Opt 

No 

No 

Surface  runofi 

UT 

No 

Opt 

No 

No 

Opt 

Seepage  fac£ 

on/off 

No 

No 

Opt 

Opt 

Opt 

40 

Hydraulic  Base  Boundary 

Constant  flux 

LH- 

Yes 

No 

Opt 

Opt 

Opt 

Constant  heac 

L 

Yes 

No 

Opt 

Opt 

Opt 

Variable  flux 

No 

No 

Opt 

Opt 

Opt 

Variable  heac 

L 

No 

No 

Opt 

Opt 

Opt 

Free  drainage 

Opt 

No 

Opt 

Opt 

No 

Deep  drainage 

No 

No 

Opt 

Opt 

No 

41 

Concentration  Surface  Boundary 

M/L3 

Inflow  solution 

M/L3 

Opt 

Opt 

Opt 

Opt 

Ambient  soli  solution 

M/L3 

Opt 

Opt 

Opt 

Opt 

42 

Concentration  Base  Bounder) 

Exit  solution 

M/L3 

Opt 

Opt 

Opt 

Opt 

Exit  flux 

M/T 

No 

Opt 

Opt 

Opt 

Convective  flow 

UT 

Opt 

No 

No 

No 

Zero  gradient 

No 

Opt 

No 

No 

Volatile  flux 

MJJ 

No 

No 

Opt 

No 

Zero  solute  flu> 

No 

No 

Opt 

Opt 

CONCENTRATION  VALUES: 

43 

Concentrations 

Location  of  maximum  concentratior 

M/L3 

Yes 

Yes 

Yes 

Yes 

Yes 

Location  of  zero  concentratior 

M/L3 

Yes 

Yes 

Yes 

Yes 

Yes 

Conceptual  initial  concentration  profile  with  deptt 

M/L3 

Yes 

Yes 

Yes 

Yes 

Yes 

Additional  2-D  parametei 

Unit 

Program 

CHEMFLO   1     SEVIEW     1  UNSATCHEM  1  HYDRUS2D  1  VS2DTI 

44 

Degree  ot  anisotropy 

No      1      Yes      1       No      |      Yes      I  Yes 

• 


Table  A2  Complete  input  data  requirements  for  CHEMFLO,  HYDRUS,  UNSATCHEM,  VS2DTI 


This  table  identifies  wtiat  parameters  need  to  be  input  with  eacli  model.  If  the  cell  says  Yes"  then  the  variable  must  be  entered  into  the 
model.  If  the  cell  says  "No"  then  there  is  no  option  to  enter  the  variable  into  the  model.  If  the  cells  sayOpt"  then  the  variable  can  be  input 
into  the  model,  but  nothing  needs  to  be  input  into  the  model  if  it  is  not  clicked  on  in  the  model. 


Must  be  measured/providec 
Must  be  estimateo 
Not  needed/Defaulte< 


Table  A1  Detailed  input  data  for  CHEMFLO,  HYDRUS,  UNSATCHEM,  VS2DTI 


VsrisbiG 

Unit 

Proc 

ram 

CHEMFLO 

HYDRUS 

UNSATCHEM 

VS2DTI 

SOIL  PROPERTIES: 

Moisture  Content 

saturated  moisture  content  (85)  -  Coarse  Material 

0.25 

0.25 

0.25 

0.25 

residual  moisture  content  (6,)  -  Coarse  Material 

0.05 

0.05 

0.05 

0.05 

saturated  moisture  content  (65)  -  Fine  Material 

0.36 

0.36 

0.36 

0.36 

residual  moisture  content  (6r)  ■  Fine  Material 

0.3 

0.3 

0.3 

0.3 

parameter  of  the  soil  retention  curve  (9p) 

No  Option 

No  Option 

(Gr) 

No  Opfion 

parameter  of  the  soil  retention  curve  (Gm) 

No  Option 

No  Option 

(Gs) 

No  Opfion 

Permeability 

saturated  hydraulic  conductivity  (Ks)  -  Coarse  Material 

m/d 

6.9 

6.9 

6.9 

6.9 

saturated  hydraulic  conductivity  (Ks)  -  Fine  Material 

m/d 

8.6x10"* 

8.6x10"^ 

8.6x10"" 

8.6x10"* 

unsaturated  hydraulic  conductivity  (Kk)  corresponding  to  Ok 

m/d 

No  Option 

No  Option 

Ks 

No  Opfion 

relative  hydraulic  conductivity  (Kr) 

m/d 

No  Option 

No  Option 

Ks 

No  Opfion 

reduction  in  hydraulic  conductivity  due  to  solution  chemistry 

on/off 

No  Option 

No  Option 

Turned  Off 

No  Opfion 

swcc 

Van  Genuchten  -  a  parameter  -  Coarse  Material 

m' 

12 

12 

12 

12 

Van  Genuchten  -  p  parameter  -  Coarse  Material 

2.7 

2.7 

2.7 

2.7 

Van  Genuchten  -  a  parameter  -  Coarse  Material 

m"' 

1.5 

1.5 

1.5 

1.5 

Van  Genuchten  -  p  parameter  -  Coarse  Material 

1.1 

1.1 

1.1 

1.1 

Pore  Conductivity  Parameter  -  Coarse  Material 

No  Option 

 ^  

No  Option 

No  Option 

Pore  Conductivity  Parameter  -  Fine  Material 

No  Option 

0.5 

No  Opfion 

No  Option 

Specific  Storage  -  Coarse  IVIaterial 

m 

No  Opfion 

No  Option 

No  Opfion 

1x10"* 

Specific  Storage  -  Fine  Material 

m'^ 

No  Option 

No  Option 

No  Opfion 

1x10-' 

Degree  of  anisotropy  -  Coarse  Material 

No  Option 

1 

No  Opfion 

1 

Degree  of  anisotropy  -  Fine  Material 

"" 

No  Option 

1 

No  Opfion 

1 

Bulk  density  -  Coarse  Material 

kg/m 

1750 

1750 

1750 

1750 

Bulk  density  -  Fine  Material 

kg/m 

1750 

1750 

1750 

1750 

SOLUTE  PROPERTIES: 

Dispersion  coefficient 

longitudinal  dispersivity 

m 

0.4 

0.4 

0.4 

0.4 

transverse  dispersivity 

m 

No  Option 

0.04 

No  Opfion 

0.04 

Diffusion  coefficient  (ions  in  soil) 

diffusion  coefficient  (Na+)  -  coarse  soil 

— m^/d 

 5  

2.9x10 

 5  

2.9x10 

 5  

2.9x10 

 E  

2.9x10 

diffusion  coefficient  (Na+)  -  fine  soil 

m  la 

No  Option 

No  Option 

4.0x10"* 

4.0x10'* 

diffusion  coefficient  (air) 

m^/d 

No  Option 

0 

Not  Required 

No  Opfion 

Kd 

uniform  soil  water  partition  coefficient 

mVkg 

0 

No  Opfion 

No  Opfion 

Turned  Off 

Freundlich  adsorption  isotherm 

m^/kg 

No  Option 

Below 

No  Opfion 

Turned  Off 

adsorption  isotherm  coefficient  kd 

m^/kg 

No  Option 

0 

No  Opfion 

No  Option 

adsorption  Isotherm  coefficient  v 

m^/kg 

No  Option 

0 

No  Option 

No  Option 

adsorption  Isotherm  exponent  p 

No  Option 

1 

No  Opfion 

No  Opfion 

Ion  Exchange 

Km:  Gaines-Thomas  ion-exchange  selectivity  coefficients  (Ca-Na) 

No  Option 

No  Option 

No 

8.5 

Km:  Gaines-Thomas  ion-exchange  selectivity  coefficients  (Ca-Mg) 

No  Option 

No  Option 

No 

1.2 

Km:  Gaines-Thomas  ion-exchange  selectivity  coefficients  (Ca-K) 

No  Option 

No  Option 

No 

0.05 

Kg:  Gapon  ion-exchange  selectivity  coefficients  (Ca-Na) 

No  Option 

No  Option 

2.9 

No  Opfion 

Kg:  Gapon  ion-exchange  selectivity  coefficients  (Ca-Mg) 

No  Option 

No  Option 

1.2 

No  Opfion 

Kg:  Gapon  ion-exchange  selectivity  coefficients  (Ca-K) 

No  Option 

No  Option 

0.2 

No  Opfion 

Cation  exchange  capacity 

meq/kg 

No  Option 

No  Opfion 

90 

90 

Calcite  surface  area 

m' 

No  Option 

No  Opfion 

0 

No  Opfion 

Dolomite  surface  area 

m' 

No  Option 

No  Opfion 

0 

No  Opfion 

Specify  kinetic  precipitation/dissolution 

on/off 

No  Option 

No  Option 

Turned  Off 

No  Opfion 

Organic  carbon  fraction 

0 

Not  Required 

Not  Required 

Not  Required 

Decay  coefficients 

first  order  decay  coefficients 

 — 

d 

0 

0 

No  Opfion 

0 

first  order  decay  coefficient  for  chain  react's 

d"^ 

No  Option 

0 

No  Opfion 

No  Option 

Production  coeffcients 

zero  order  production  coefficient(gas) 

 ■  

d' 

No  Option 

0 

Turned  off 

No  Option 

zero  order  production  coefficient  (solid) 

d 

0 

0 

No  Otpion 

No  Option 

Chemical  Non-Equilibrium  Parameters 

Fract 

No  Option 

1 

No  Opfion 

No  Option 

Alpha 

d-^ 

No  Option 

0 

No  Opfion 

No  Option 

Physical  Non-Equilibrium  Parameters 

Thimob 

m^/m^ 

No  Option 

0 

No  Opfion 

No  Option 

Temperature  Dependent  Parameters 

Henry"' 

No  Option 

0 

No  Opfion 

Not  Required 

temperature  dependence  for  water  flow  parmeters 

on/off 

No  Option 

Turned  Off 

No  Opfion 

No  Opfion 

(1)  Simunek,  Suarez,  and  Sejna  (1996)  suggest  that  the  hydraulic  characterisfics  contain  9  unknown  parameters,  B„  Gs,  Gp,  G^,  a,  p,  Kj,  Kk,  G^ 
When  %=6„  Gm=G|<=Gs,  and  K,,  =  Kg,  the  soil  hyddraulic  functions  reduce  to  the  original  expressions  of  van  Genuchten  (1980). 

(2)  Simunek,  Sejna  and  Van  Genuchten  (1999)  suggest  that  the  pore-connectivity  parameter  (I)  in  the  hydraulic  conducfivity  function  was 
esfimated  (Mualem,  1976)  to  be  about  0.5  as  an  average  for  most  soils. 

(3)  Fract.  -  Dimensionless  fracfion  [-]  of  the  sorpfion  sites  classified  as  type-1,  (i.e.,  sites  subject  to  instantaneous  sorpfion)  when  chemical 
nonequilibrium  is  simulated,  or  dimensionless  fracfion  of  sorpfion  sites  in  contact  with  mobile  water  when  physical  nonequilibrium  is  simulated. 
Author  says  default  is  one  if  non-equilibrium  is  not  considered. 

(4)  Alpha  -  First-order  rate  transfer  coefficient  for  nonequilibrium  sorption  when  chemical  nonequilibrium  is  simulated,  or  for  exchange  between 
mobile  and  immobile  liquid  regions  when  physical  nonequilibrium  is  simulated 

(5)  Thimob  -  Immobile  water  content  when  physical  nonequilibrium  is  simulated. 

(6)  Henry  -  Equilibrium  distribufion  constant  between  liquid  and  gas  phases 
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Figure  B1  -  Source  applied  to  a  coarse-grained  soil  column  for  5  Years. 


Figure  B2  -  Source  applied  for  to  a  coarse-grained  soil  column  10  years. 


CHEMFLO  -^VS2DTI  ^HYDRUS  -b-UNSATCHEM 


Figure  B3  -  Source  applied  to  a  coarse-grained  soil  column  for  15  years. 


10  

? 

£  15  

o 

a 

20  

25  

30  I  \  \  \  ^  1  1  \  \  \  

OT-cN|cO'si;iocqtv.c3qa> 
oooodoood 

C/Co 

-^CHEMFLO  -^VS2DTI  -^HYDRUS  -b-UNSATCHEM 


Figure  B-4  Source  applied  to  a  coarse-grained  soil  column  for  20  years. 


Figure  B5  -  Source  removed  from  a  coarse-grained  soil  column  for  5  years. 


Figure  B6  -  Source  removed  from  a  coarse-grained  soil  column  for  10  years. 
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Figure  B7  -  Source  removed  from  a  coarse-grained  soil  column  for  15  years. 


Figure  B8  -  Source  removed  from  a  coarse-grained  soil  column  for  20  years. 


Figure  B9  -  Source  added  to  a  fine-grained  soil  column  for  5  years. 
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Figure  B10  -  Source  added  to  a  fine-grained  soil  column  for  10  years. 
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Figure  B11  -  Source  added  to  a  fine-grained  soil  column  for  15  years. 


C/C, 


CHEMFLO  -^VS2DTI  -a-HYDRUS  -b-UNSATCHEM 


Figure  B12  -  Source  added  to  a  fine-grained  soil  column  for  20  years. 
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Figure  B13  -  Source  removed  from  a  fine-grained  soil  column  for  5  years. 
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Figure  B14  -  Source  removed  from  a  fine-grained  soil  column  for  10  years. 
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Figure  B15  -  Source  removed  from  a  fine-grained  soil  column  for  15  years. 
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Figure  B16  -  Source  removed  from  a  fine-grained  soil  column  for  20  years. 
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Figure  B17  -  Source  added  to  a  fine-coarse-fine  layered  column  for  5  years. 
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Figure  B18  -  Source  added  to  a  fine-coarse-fine  layered  column  for  10  years. 
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Figure  B19  -  Source  added  to  a  fine-coarse-fine  layered  column  for  15  years. 
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Figure  B20  -  Source  added  to  a  fine-coarse-fine  layered  column  for  20  years. 
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Figure  B21  -  Source  removed  from  a  fine-coarse-fine  layered  column  for  5  years 
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Figure  B22  -  Source  removed  from  a  fine-coarse-fine  layered  column  for  10  years. 
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Figure  B23  -  Source  removed  from  a  fine-coarse-fine  layered  column  for  15  years. 
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Figure  B24  -  Source  removed  from  a  fine-coarse-fine  layered  column  for  20  years. 
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Figure  B25  -  Source  added  to  a  coarse-fine-coarse  layered  column  for  5  years. 
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Figure  B26  -  Source  added  to  a  coarse-fine-coarse  layered  column  for  10  years. 
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Figure  B27 


Source  added  to  a  coarse-fine-coarse  layered  column  for  15  years. 


Figure  B28  -  Source  added  to  a  coarse-fine-coarse  layered  column  for  20  years. 
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Figure  B29  -  Source  removed  from  a  coarse-fine-coarse  layered  column  for  5 

years. 


Figure  B30  -  Source  removed  from  a  coarse-fine-coarse  layered  column  for  10 

years. 


Figure  B31  -  Source  removed  from  a  coarse-fine-coarse  layered  column  for  15 

years. 


Figure  B32  -  Source  removed  from  a  coarse-fine-coarse  layered  column  for  20 

years. 
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Figure  B34  -  Distributed  source  on  a  coarse-fine-coarse  layered  column  year  5. 
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Figure  B35  -  Distributed  source  on  a  coarse-fine-coarse  layered  column  year  10. 


Figure  B36  -  Distributed  source  on  a  coarse-fine-coarse  layered  column  year  20. 
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This  information  was  prepared  by  MDH  Engineered  Solutions  Corp.  (MDH)  to  provide  a 
qualitative  comparison  of  the  programs  HYDRUS  and  VS2DTI.  Each  program  outputs 
the  concentration  data  in  a  different  format.  In  each  program,  red  indicates  a  C/Co  of 
one  and  blue  indicates  a  C/Co  of  zero.  However,  the  variation  in  contour  colours 
between  C/Co  =  1  and  C/Co  =  0  are  not  the  same.  The  information  provided  in  this 
appendix  can  be  used  to  compare  the  extent  of  the  plume  and  the  position  of  the  peak 
concentration  at  different  times. 

The  output  for  the  pipeline  scenario  show  that  both  codes  produced  comparable  results. 
The  extent  of  the  plume  and  the  location  of  the  maximum  concentration  appear  to 
broadly  the  same  for  both  programs. 

The  output  from  the  salt  pile  scenario  also  show  that  both  codes  produced  comparable 
results.  The  extent  of  the  plume  and  the  location  of  the  peak  concentration  appear  to 
broadly  similar  after  20  years  of  source  application. 


Figure  C1  -  2D  plume  for  pipe  line  after  1  year  of  source  application  using  IHYDRUS. 


^  iTSm  ^ 

Figure  C2  -  2D  plume  for  pipe  line  after  1  year  of  source  application  using  VS2DTI. 


Figure  C3  -  2D  plume  for  pipe  line  after  source  removed  for  5  years  using  HYDRUS. 


^  11.5m  ^ 

Figure  C4  -  2D  plume  for  pipe  line  after  source  removed  for  5  years  using  VS2DIT. 


Figure  C7  -  2D  plume  for  salt  pile  after  20  years  of  source  application  using  HYDRUS. 


Figure  C11  -  2D  plume  for  salt  pile  after  source  removed  for  20  years  using  HYDRUS. 
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The  output  from  the  salt  pile  scenario  completed  using  VS2DTI  provided  in  this  appendix 
show  the  results  of  a  case  when  numerical  "artifacts"  were  detected  in  a  simulation. 
When  compared  with  the  VS2DTI  results  from  Appendix  C,  the  extent  of  the  plume  and 
the  location  of  the  peak  concentration  appear  to  broadly  similar  after  20  years  of  source 
application.  However,  the  results  provided  in  this  appendix  illustrate  that  after  source 
cessation  the  simulation  does  not  match  the  solution  provided  in  Appendix  C.  The 
reason  for  the  mismatch  is  that  the  VS2DTI  results  are  incorrect  because  the  grid 
spacing  and  time  stepping  did  not  satisfy  necessary  Peclet,  Courant  and  Von  Neumann 
conditions  to  suppress  numerical  "artifacts",  including  numerical  oscillation  and 
numerical  dispersion,  for  the  later  part  of  the  solution  process.  This  illustrates  the 
difficulties  that  might  be  encountered  if  modelling  is  carried  out  by  inexperienced 
practitioners. 


Figure  D3  -2D  plume  for  salt  pile  after  source  removed  for  20  years. 
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